{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# UNDERSTANDING EXPERIMENTAL DATA"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This chapter is about `understanding experimental data`. \n",
    "\n",
    "We will\n",
    "\n",
    "* make extensive use of `plotting to visualize the data`, and \n",
    "\n",
    "* show how to use `linear regression` to build a model of experimental data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## The Behavior of Projectiles\n",
    "\n",
    "We decided to use one of our springs to build `a device capable of launching a projectile`. \n",
    "\n",
    "We used the device **four** times to fire a projectile at **a target** 1080 inches from the **launching point**. \n",
    "\n",
    "Each time, we measured `the height of the projectile` at various `distances from the launching point`.\n",
    "\n",
    "* The `launching point` and the `target` were at `the same height`, which we treated as `0.0` in our measurements.\n",
    "\n",
    "![projectile](./img/projectile.jpg)\n",
    "\n",
    "The data was stored in a file `./data/projectileData.txt`\n",
    "\n",
    "* The `first column` contains `distances of the projectile from the target`. \n",
    "\n",
    "* The `other columns` contain `the height of the projectile at that distance` for each of the four trials. \n",
    "\n",
    "All of the measurements are in **inches**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting ./data/projectileData.txt\n"
     ]
    }
   ],
   "source": [
    "%%file ./data/projectileData.txt\n",
    "Distance  \ttrial1\ttrial2\ttrial3\ttrial4\n",
    "1080\t  \t0.0\t    0.0\t\t0.0\t\t0.0\n",
    "1044\t\t2.25\t3.25\t4.5\t\t6.5\t\n",
    "1008\t\t5.25\t6.5\t\t6.5\t\t8.75\n",
    "972\t\t\t7.5\t\t7.75\t8.25\t9.25\n",
    "936\t\t\t8.75\t9.25\t9.5\t\t10.5\n",
    "900\t\t\t12.0\t12.25\t12.5\t14.75\n",
    "864\t\t\t13.75\t16.0\t16.6\t16.5\n",
    "828\t\t\t14.75\t15.25\t15.5\t17.5\n",
    "792\t\t\t15.5\t16.0\t16.6\t16.75\n",
    "756\t\t\t17.0\t17.0\t17.5\t19.25\n",
    "720\t\t\t17.5\t18.5\t18.5\t19.0\n",
    "540\t\t\t19.5\t18.5\t19.0\t19.0\n",
    "360\t\t\t18.5\t18.5\t19.0\t19.0\n",
    "180\t\t\t13.0\t13.0\t13.0\t13.0\n",
    "0\t\t\t0.0\t\t0.0\t\t0.0\t\t0.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following code was used to plot :\n",
    "\n",
    "* the **mean altitude of the projectile(Y)** against the **distance from the point of launch(X)**.\n",
    "\n",
    "* the best `linear` and `quadratic` fits to the points.\n",
    "\n",
    "```python\n",
    "d2h={'d':None,'h':[]}\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "def getTrajectoryData(fileName):\n",
    "    dataFile = open(fileName, 'r')\n",
    "    d2hs=[]\n",
    "    discardHeader = dataFile.readline()\n",
    "    for line in dataFile:\n",
    "        dhs = line.split()\n",
    "        d2h={'d':None,'h':[]}\n",
    "        # the distance in first column\n",
    "        d2h['d']=float(dhs[0])\n",
    "        trials=len(dhs)-1\n",
    "        for i in range(trials):\n",
    "            d2h['h'].append(float(dhs[i+1]))\n",
    "        d2hs.append(d2h)\n",
    "    dataFile.close()\n",
    "    return  d2hs,trials"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'d': 1080.0, 'h': [0.0, 0.0, 0.0, 0.0]}\n",
      "{'d': 1044.0, 'h': [2.25, 3.25, 4.5, 6.5]}\n",
      "{'d': 1008.0, 'h': [5.25, 6.5, 6.5, 8.75]}\n",
      "{'d': 972.0, 'h': [7.5, 7.75, 8.25, 9.25]}\n",
      "{'d': 936.0, 'h': [8.75, 9.25, 9.5, 10.5]}\n",
      "{'d': 900.0, 'h': [12.0, 12.25, 12.5, 14.75]}\n",
      "{'d': 864.0, 'h': [13.75, 16.0, 16.6, 16.5]}\n",
      "{'d': 828.0, 'h': [14.75, 15.25, 15.5, 17.5]}\n",
      "{'d': 792.0, 'h': [15.5, 16.0, 16.6, 16.75]}\n",
      "{'d': 756.0, 'h': [17.0, 17.0, 17.5, 19.25]}\n",
      "{'d': 720.0, 'h': [17.5, 18.5, 18.5, 19.0]}\n",
      "{'d': 540.0, 'h': [19.5, 18.5, 19.0, 19.0]}\n",
      "{'d': 360.0, 'h': [18.5, 18.5, 19.0, 19.0]}\n",
      "{'d': 180.0, 'h': [13.0, 13.0, 13.0, 13.0]}\n",
      "{'d': 0.0, 'h': [0.0, 0.0, 0.0, 0.0]}\n"
     ]
    }
   ],
   "source": [
    "fileName='./data/projectileData.txt'\n",
    "d2hs,trials = getTrajectoryData(fileName)\n",
    "for item in d2hs:\n",
    "    print(item) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1 Linear Regression to Find a Fit\n",
    "\n",
    "Whenever we fit any curve (including a line) to data we need some way to decide `which curve is the best fit for the data`. \n",
    "\n",
    "This means that we need to define <b style=\"color:blue\">an objective function</b> that provides `a quantitative assessment of how well the curve fits the data`. \n",
    "\n",
    "Once we have such `a function`, finding the best fit can be formulated as\n",
    "\n",
    "* finding a curve that <b style=\"color:blue\">minimizes (or maximizes)</b> the value of that function, i.e., as **an optimization problem** \n",
    "\n",
    "The most commonly used objective function is called <b style=\"color:blue\">least squares(最小二乘)</b>,\n",
    "\n",
    "The objective function is then defined as\n",
    "\n",
    "$$\\sum_{i=0}^{len(observed)-1}(observed[i]-predicted[i])^2$$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Numpy.polyfit**\n",
    "\n",
    "http://docs.scipy.org/doc/numpy/reference/generated/numpy.polyfit.html\n",
    "\n",
    "`Numpy` provides a function, `polyfit`, that finds the best Least squares polynomial fit.\n",
    "```python\n",
    "numpy.polyfit(x, y, deg, rcond=None, full=False, w=None, cov=False)\n",
    "```\n",
    "The algorithm used by `polyfit` is called <b style=\"color:blue\">linear regression</b>.\n",
    "\n",
    "Fit a polynomial of degree `deg` to points (x, y). \n",
    "\n",
    "$$p(x) = p[0] * x^{deg} + ... + p[deg]$$\n",
    "\n",
    "Returns a vector of **coefficients `p`** that minimises the squared error.\n",
    "\n",
    "\n",
    "\n",
    "```python\n",
    "numpy.polyfit(observedXVals, observedYVals, n)\n",
    "numpy.polyfit(observedXVals, observedYVals, 1) # y = ax + b\n",
    "numpy.polyfit(observedXVals, observedYVals, 2) # y = ax^2 + bx+c\n",
    "```\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0.16071429,  0.50071429,  0.22142857])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "x = np.array([0.0, 1.0, 2.0, 3.0,  4.0,  5.0])\n",
    "y = np.array([0.0, 0.8, 0.9, 0.1, -0.8, -1.0])\n",
    "z = np.polyfit(x, y, 2)\n",
    "z"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2 Using numpy.polyfit to fit the experimental data\n",
    "\n",
    "* linear fit : `a,b = np.polyfit(distances, meanHeights, 1)`\n",
    "  \n",
    "* quadratic fit: `a,b,c = np.polyfit(distances, meanHeights, 2)`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def processTrajectories(d2hs,trials):\n",
    "    numTrials = trials\n",
    "    distances=[]\n",
    "    meanHeights=[]\n",
    "    \n",
    "    for item in d2hs:\n",
    "        distances.append(item['d'])\n",
    "        meanHeights.append(np.mean(item['h']))\n",
    "    \n",
    "    plt.title('Trajectory of Projectile (Mean Height of {} Trials'.format(trials))\n",
    "    plt.xlabel('Inches from Launch Point')\n",
    "    plt.ylabel('Mean Height of the projectile')\n",
    "    \n",
    "    # the experimental data\n",
    "    plt.plot( distances, meanHeights, 'ro')\n",
    "    \n",
    "    # Linear Fit\n",
    "    a,b = np.polyfit(distances, meanHeights, 1)\n",
    "    altitudes = a*np.array(distances) + b\n",
    "    plt.plot(distances, altitudes, 'g', label = 'Linear Fit')\n",
    "    \n",
    "    # Quadratic Fit\n",
    "    a,b,c = np.polyfit(distances, meanHeights, 2)\n",
    " \n",
    "    altitudes = a*(np.array(distances)**2) +  b*np.array(distances) + c\n",
    "    plt.plot(distances, altitudes, 'b:', label = 'Quadratic Fit')\n",
    "    plt.legend(loc=\"best\")\n",
    "    plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEWCAYAAABhffzLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAABB6klEQVR4nO3deXwU9f348dc7IRAS7lNUSBARRBBEQCxa8b6lVesVtZ60Wms9aj3wp/bA6ldaPGul1uIRb8W7qCCo9QQUEUEEMRwKCSJyJFxJ3r8/PrNh2exuNsnuzh7v5+Oxj+wcO/Oenc1nZj7zmfdHVBVjjDHZI8fvAIwxxiSXFfzGGJNlrOA3xpgsYwW/McZkGSv4jTEmy1jBb4wxWcYK/gQSkf+KyC/9jiOeROQSESkXkU0i0jmB6zlYRBYlerkiUiYiRzRxWa1EZIGI9IhfhKlFRP4pIv8vxnkni8hfEhjLz0Vkhffb2y9R6wmz3k0iskcM8xWLiIpIi2TE1RxW8IfwdnLgVSsim4OGSxqzLFU9VlUfbmY8t4jIY81ZRryISB7wd+AoVW2jqmtDpgd++IHvq0xErmvKulT1XVXtF4eYVUT2jPdyPWOBd1R1lbeuyd76xoTEMNEbf16c1huz0O33xsX8m1LVX6vqnxMVSyNNAC7zfnufRlnPId66wh6EvBOywG90u4hsCxr+Z+j83vqWNiPulJPyR6ZkU9U2gfciUgZcpKrTQucTkRaqWp3M2JoiznF2B/KBLxqYr4OqVovIgcB0EZmrqlMTGJdffo0r/IN9BZwLvAhuO4HTgK+TG1pGKqKB3553cnIX8FGkeVT12KD5JwMrVfXGMMvKhN9oWHbGHyMRGS0iK0XkWhFZDfxHRDqKyCsiskZE1nnvdw/6zEwRuSho+AIRWejN+7qIFAVN20dE3hSRH7yqlBtE5BjgBuB072zkM2/eXUXkJW/eJSJycdBybhGRZ0XkMRHZAFwnIlXB1TIiMtSLOS/MdrYSkTtF5Dvvdac3bi8gUEXyo4i81dB3pqof4P5RB0b4/sKuK/j7DoprVxF5zov7GxG5PGharvd9fS0iG0Vkjoj0FJF3vFk+876/00OXG7LtOSJynbectSLytIh0ijBvL2AP6hcwLwMHiUhHb/gYYB6wOuTz0X4Ld4mr0tjgbcvBQdNu8eJ6xNvWL0RkWJTd0CAR6R/021skIqcFTdup+kZE/iAiq7z9dZHUP4vvKCKverF9JCJ9vM/V2xdh4sgRkRtFZJmIVHjb2N77nWwCcr3PRzuIXg28AXzZxO9CReQ3IrIYWBw0bk/v/fEi8qm3b1aIyC1RlnWeiCz1votvpJE1BolkBX/j7AJ0wp15jMV9f//xhnsBm4F7w31Q3OX/DcDJQFfgXeAJb1pbYBowFdgV2BOY7p0l3wo85V1uDvYW9ySw0pv3VOBWETksaHVjgGeBDsDfgJm4s86Ac4AnVXV7mFDHASOBIcBgYARwo6p+BezjzdNBVQ8L89ng7RURGeV9JnBZHvr9hV1XmGXl4ArUz4DdgMOBK0TkaG+Wq4AzgeOAdsAFQJWq/tSbPtj7/p6KFjPwW+BnwCG473YdcF+EeQcBS8OcEW7Bne2f4Q2fCzwSsj0RfwueWbjvpBPwOPCMiOQHTT8J9xvoALxEhN9cLESkEHjTW083L+5/iMiAMPMeg/uuj8D9RkeHWeQZwB+BjsASYDxAjPviPO91KO6g2ga4V1W3Bl2JD1bVPhG2pQi37/8UdaMb9jPgAKDedwBU4vZpB+B44BIR+VmYWAqBu4FjVbUt8BNgbjPjih9VtVeEF1AGHOG9Hw1sA/KjzD8EWBc0PBNXVQTwX+DCoGk5QBWuEDwT+DTCMm8BHgsa7gnUAG2Dxv0VmBw0/zshyzgdeM97n4s7+xwRYX1fA8cFDR8NlHnviwEFWkT4bGD6j7hCcyFweaTvr4F1jcZdgoP7J1wesq7rgf947xcBYyLEpMCeQcN1yw2zjxcChwdN6wFsD7e9QAnwYci4ycBfgIOAD3CFQznQGvgfcF5Dv4UI27AOV+AF9u+0oGkDgM1RfpMKbPD2SeC1JfCb8n4b74Z85gHg5uBt8t4/BPw1aL49g79fb94Hg6YfB3wZaV+EiXU6cGnQcL/g7z+Gz78InB4ad7RX6HzeOg6L9hsKmXYnMDH0/wMo9L7rU4DWDcWR7Jed8TfOGlXdEhgQkQIRecC7NN0AvAN0EJHcMJ8tAu4SkR9F5EfgB0BwZ7A9ib0OeFfgB1XdGDRumbecgBUhn3kRGCAivYEjgfWq+nGU5S8LWfauMcYW0EVVO6rq3qp6d9D4nb6/RqyrCNg18N15398NuHsO0LjvL5oiYErQOhbiDrLdw8y7DmgbbiGq+j/cmfw44BVV3RxmPZF+C4jI771qoPXe9PZAl6DPB1cbVQH5Er0lyVBV7RB4AbeFxHJAyHdbgrs6C7UrO/+2Qn9n4WJrE2aeSML9HloQ/vvfiYiciDsZauiqLhbhtiuwngNEZIa4Ksf1uPs8XULnU9VK3EH118Aqr/qrfxxiiwu7uds4oalMr8adlRygqqtFZAiuWkPCfHYFMF5VS0MneJeoZ9T/SNh1fgd0EpG2QYV/L+DbSJ9R1S0i8jRwNtAfeDTCugLLD76J1ssbFw/htiWWda0AvlHVvhGWuwLoA8xvZnwrgAtU9b0Y5p0H9JbINwAfA27CVVuEW0+k38LBwB9w1VlfqGqtiKwj/G8qHlYAb6vqkTHMuwrYPWi4Z5xjCfweAnoB1birpoYcDgwTd/8I3MGyRkQGqeqYKJ8LJ1rK4sdxVWvHev9XdxKm4AdQ1deB10WkNe5K8F/AweHmTTY742+etrh6/R/F3QS8Ocq8/wSuF5F9ALybVr/wpr0C9BCRK7wbWW1F5ABvWjlQ7NVzo6orgPeBv4pIvojsC1yIK2iieQRXf3oS0Qv+J4AbRaSriHTBFV6Jak4a67o+BjaKuzHcWtzN3IEiMtyb/iDwZxHp691b2Fd23Mwux9UXx+KfwHjvQIwXV9hCQ1VX4uqwR0RY1t24q6t3wkyL9ltoiyvs1gAtROQm3H2LRHkF2EtEzhGRPO81XET2DjPv08D5IrK3iBQAMbXvD9LQvngCuFJEeotIG3bc34qlZc3/A/bCVbcOwd37+BdwfiNjbEhb3BX3FhEZAZwVbiYR6S4iY7y6/q3AJqA2zrE0mRX8zXMnrv72e+BD3M3ZsFR1CnA78KRXLTQfONabthFXSJyIu1RezI4zxWe8v2tF5BPv/Zm4+sTvgCm4+th6TU5D1v8e7of3iaouizLrX4DZuDPaz4FPvHGJENO6VLUGOAH3D/0N7vt+EHdWB+7ZgqdxrTk2AP/G7RdwdeIPe9UYwTe4w7kLV2C8ISIbcfv0gCjzP4C7UV6Pqv6gqtPVq/wNmRbxtwC8jvsdfYWr6thClKqH5vJ+e0fhrji/w/3+bgdahZn3v7gD2gzcQe9Db9LWGFd3C9H3xUO4k5J3cPt5C+6Ge0zboaqrAy/cCVmlqv4QY2yxuhT4k/f7uAn3uwsnB3cj/DtcVd4hwCVxjqXJJMzv0sSJuCZsD6rqIw3OnATimmA+rqoP+h1LQ7xWSg+qaqxn60knrunpp7gbwqv8jifZvKuC+UCrGM/KTYqwM/4E8S6F98CdufjOqxYZCsTj5lcyDCRFvrtI1DUzHJBNhb64tAmtxD2ncDvwshX66ccK/gQQkW64S+a3cc34fCUiD+OeE7gipDVQShKRu4Arce3BTWr5FVCBa0VVQwpVX5jYWVWPMcZkGTvjN8aYLJMW7fi7dOmixcXFfodhjDFpZc6cOd+ratfQ8WlR8BcXFzN79my/wzDGmLQiImGbbltVjzHGZBkr+I0xJstYwW+MMVkmLer4jTGpZ/v27axcuZItW7Y0PLNJqPz8fHbffXfy8ur1rRSWFfzGmCZZuXIlbdu2pbi4GJFEJQ81DVFV1q5dy8qVK+ndu3dMn7GqHmNMk2zZsoXOnTtboe8zEaFz586NuvKygt8Y02RW6KeGxu4HK/iNaaLZs+Gee2DbNr8jMaZxrOA3phFqg7rSePBBuPdeCNxPW7AAtofrvt4kTJs29Xt2/Oc//8kjjyQ3E/ro0aPp168fQ4YMYciQITz77LNcdNFFLFiwAIBbb701qfE0JC2StA0bNkztyV3jt/feg7POgjffhL32gooKEIGuXaG6Gnr1gkMPhdJ6HSpmpoULF7L33uE66kqeNm3asGnTpqSus67D8pwd582jR49mwoQJDBs2LOxnkhFnuP0hInNUtV5QdsZvTBRffAELF7r3e+4Je+8NVVVuuFs3V+iDOwA8+CD81usv6ocf4Kij4ONIXdqbhLnllluYMGEC4Arka6+9lhEjRrDXXnvx7rvvAlBTU8M111zD8OHD2XfffXnggQcA2LRpE4cffjhDhw5l0KBBvPjiiwCUlZXRr18/zj33XAYOHMiKFQ13ijZ69Ghmz57Nddddx+bNmxkyZAglJSUJ2urGseacxkSwbRsccggcdhg8/TR07w5TI3SumZsLxx23Y3jpUvj6a2jZ0g2vWAHr18PAgYmP2w9XTL2CuavnxnWZQ3YZwp3H3Nns5VRXV/Pxxx/z2muv8cc//pFp06bx73//m/bt2zNr1iy2bt3KqFGjOOqoo+jZsydTpkyhXbt2fP/994wcOZKTTjoJgMWLF/Pwww8zcuTIsOspKSmhdWvX4+f06dPrxt92223ce++9zJ07t9nbEi9W8BsT5Mkn4ZVX4LHHXKH97LNNK6yHDYMlS9yVAMBdd7lXeTl06hTfmE10J598MgD7778/ZWVlALzxxhvMmzePZ599FoD169ezePFidt99d2644QbeeecdcnJy+PbbbykvLwegqKgoYqEPUFpaGrGqJ9UkrOAXkZ7AI0B3QIFJqnqXiHTCdf9XDJQBp6nqukTFYUxDysqgZ0931r5mDSxe7M7O27eH0aObvtzgFnbXXQc//emOQv/886FHD0ixe35NFo8z80Rp1cr1G5+bm0t1teslUlW55557OProo3ead/LkyaxZs4Y5c+aQl5dHcXFxXfv4wsLC5AaeQIms468GrlbVAcBI4DciMgC4Dpiuqn2B6d6wMb54/33o0wdeftkNX3opfPSRK/TjqUsX8GoMUHUtgXJzd0x/8kl3sDHJcfTRR3P//fez3WuG9dVXX1FZWcn69evp1q0beXl5zJgxg2XLwmY1brS8vLy6daWChJ3xex1Qr/LebxSRhcBuwBhgtDfbw8BM4NpExWFMsNpaePRRaNsWTj4ZRoyAP/0Jhg9304ML40QRgUmTdgwvWABnngl33+1uDtfWunns2aiGVVVVsfvuu9cNX3XVVTF97qKLLqKsrIyhQ4eiqnTt2pUXXniBkpISTjzxRAYNGsSwYcPo379/XOIcO3Ys++67L0OHDqU0BZp9JaU5p4gUA+8AA4HlqtrBGy/AusBwyGfGAmMBevXqtX+8jrwmO23dCq1aubPtAw6AXXeFF17wO6od5sxxVx4dOsCLL8L118Orr0KMqVd8kQrNOc0OKdWcU0TaAM8BV6jqhuBp6o46YY88qjpJVYep6rCuXev1HGZMzO65xxWqW7a4s+hXXoEpU/yOamf77+8KfYDCQujb1913ANeSaNo0d9AyJh4SWvCLSB6u0C9V1ee90eUi0sOb3gOoSGQMJvtUV8Nzz7kbtQCDB7tqneD296lcjXLEEe6sv4VXEXvrrXDDDTti3rzZv9hMZkhYwe9V4/wbWKiqfw+a9BLwS+/9L4EXExWDyU5ffw2nnuqaZIJrTXP33enbjPKNN+CJJ9z7rVtd9c8dd/gbk0lviTzjHwWcAxwmInO913HAbcCRIrIYOMIbNqZZrr4arrnGve/XD959Fy6/3N+Y4iU/31VVgauuOu88d1Ma3HMBN93k/hoTq0S26vkfEOmC+vBErddkB1X49FMYOtQNV1XtqBoBOOggf+JKtPbt4bagU6WZM2H8eJdDqHt3WLfO3SMIPDFsTDiWq8ekpb/9zT0du3SpG/7HP9xN3Gxz+unw7bcQaHU4bpxLIJdCTcZNCrKC36SFdevg2mth1iw3fNZZ8MgjEGjCnfSbtaWlUFwMOTnur49ts3fZZcf7k0+Gq67akSr6xhvdjeJMtXLlSsaMGUPfvn3ZY489uOyyy9i6dWtclj1z5kxOOOGERn2mrKyMxx9/vG549uzZXN6IOsekpXcOpBhN5df++++vJvvU1qr++KN7v3GjaufOqnfe6W9Mqqr62GOqBQWqrsbJvQoK3PgUsnmzat++quPGueHaWtWFC+O3/AULFsRvYU1QW1urw4cP14ceekhVVaurq/WCCy7Qyy+/PC7LnzFjhh5//PH1xm/fvr3Rn4nVIYccorNmzYo4vbCwMOK0cPsDmK1hylQ74zfJ1Ygz5dNO25HmoE0bl1Pnd79LRpANGDduR9vQgKoqNz6F5OfDl1/uCOujj1xa6eee8zeueHnrrbfIz8/n/PPPB1wunokTJ/LII4+wadMmJk+ezGWXXVY3/wknnMDMmTMBuOSSSxg2bBj77LMPN998c908U6dOpX///gwdOpTnn3++bvwtt9zCOeecw6hRozjnnHMoKyvj4IMPZujQoQwdOpT3338fgOuuu453332XIUOGMHHixJ2uGjZt2sT555/PoEGD2HfffXkuxh2RiPTOVvCb5CkthbFjYdkyd568bJkb9gr/igqXwTLQy9WYMXDGGTseXArT2ZI/li9v3Hgf5eSAlymYvfaCv//d9RMA7snlX/0KNmyI+PFGGT0aJk9277dvd8OBJrVVVW74qafc8Pr1bjhQtn7/vRsO5Exavbrh9X3xxRfsv//+O41r164dxcXFLFmyJOpnx48fz+zZs5k3bx5vv/028+bNY8uWLVx88cW8/PLLzJkzh9UhQSxYsIBp06bxxBNP0K1bN958800++eQTnnrqqbrqnNtuu42DDz6YuXPncuWVV+70+T//+c+0b9+ezz//nHnz5nHYYYeFja2kpKSuqmft2rV142+77TZat27N3Llzm532wdIym+SJcKasN4xDSkqYPh2uuAJGjnRpFc4+25coG9arlztohRufwjp1guCy6OuvXa9igQPqhg3Qrp0/sSXb008/zaRJk6iurmbVqlUsWLCA2tpaevfuTd++fQE4++yzmRSUVOmkk06qy7e/fft2LrvsMubOnUtubi5fffVVg+ucNm0aTz75ZN1wx44dw86XjPTOVvCb5Ak5I15HB07meUqWP85FuIeu9ttvRwuVlDV+vLtSCT6IFRS48Wnk6qvdgTYnx51xH3CAuwL4wx+atjyvFgVwN5eDhwsKdh5u337n4S5ddh4OvmEdyYABA+ry6Qds2LCB1atX069fP+bPn09tUCfJgfTK33zzDRMmTGDWrFl07NiR8847r25aNMFpmSdOnEj37t357LPPqK2tJT8/v+GAU4hV9Zjk6dULBebgGt934EfasYH8zu4fKi8vDQp9gJISl16zqMg1JyoqcsMp0q1eY+TmAqWldNi/D8ctvZfRE8ekTafBhx9+OFVVVXUdq9fU1HD11Vdz2WWX0bp1a4qLi5k7dy61tbWsWLGCj71+MDds2EBhYSHt27envLyc//73vwD079+fsrIyvv76awCeCDwuHcb69evp0aMHOTk5PProo9TU1ADQtm1bNm7cGPYzRx55JPfdd1/d8Lp1je+GJF7pna3gN8kzfjwT8m7geF5lKy0R4MWCszj7ruF+R9Z4JSXubnNtrfubhoU+UHffpcXypdzDbxmx+iUYO5ZXfj+TJPdh3mgiwpQpU3j22Wfp27cvnTt3Jicnh3He3exRo0bRu3dvBgwYwOWXX85Q72m/wYMHs99++9G/f3/OOussRo0aBUB+fj6TJk3i+OOPZ+jQoXTr1i3iui+99FIefvhhBg8ezJdffll3NbDvvvuSm5vL4MGDmThx4k6fufHGG1m3bh0DBw5k8ODBzJgxo9HbHEjv3Oy+e8M19Um1lzXnzByL7nhRb2h3j9YiqkVFKdcEMusUFe3cLBV0ObtrHlv1mmuif9Tv5pyh3nvvPe3Vq5fOmTPH71B80ZjmnEnJx99cw4YN09mzZ/sdhmmGzZt3tC4xKSQnJ2y+5xkcyohNbxGtt0HLx59aUiofvzFVVfCTn8Att/gdiaknQkukQ4uWUlgI27a53sG86nGTIazgNwmXl+dSIwcySpoUMn68a3ITLKiFUnm5K/QXLQr/8XSoMcgGjd0PVvCbhKqtdQX/XXfBccf5HY2pp4EWSj17wvz5cM45QGkpFT33r3vqOn/tWtauXWuFv89UlbVr1zaqSanV8ZuEefdd98DQ88+n/LNNpiGlpSy+6HZGbHmbv3I9v+YBtu+6KyufeootnTv7HV3Wy8/PZ/fddycvkJ3PE6mOv8EHuESkALga6KWqF4tIX6Cfqr4Sr6BNZtq2zeWFb9/e70hMs40bR/GWbzmf/3AcrwGQ99139D77bNec1aSVBs/4ReQpYA5wrqoO9A4E76vqkCTEB9gZfzpTTe3+bU2MQlr/uAfx9meYfLIjuZJJOc1p1dNHVf8P2A6gqlVE7lnLGK68ckeyLiv0M0RIXd2jnMNwZvNOt1N9Csg0Ryy5eraJSGvcQR4R6QPEp6cDk3G2bIF586zrv4wTkp/oNJ5mY15nDrpjjM+BmaaIpeC/GZgK9BSRUlwn6uclMiiTvvLz4fXX/Y7CxF0gRcC4cbB8Ofm9duE344dBSQnffw8LFrgmuyY9xNSqR0Q6AyNxVTwfqur3iQ4smNXxp741a+Dmm+H226FtW7+jMclUUgJTp8I332RPWud00eg6fhEZGngBRcAq4DuglzfOmDpvv+36wPUSG5oscuedrl9fK/TTR7Sqnr9FmaZA+O5jTFY69VQ45BDo2tXvSEyyde26Y7//97/uoa+BA/2NyUQXseBX1UOTGYhJTy+84DrNGDnSCv1st3UrXHIJDBq0owtFk5oiFvwicpiqviUiJ4ebrqrPhxtvskdNDdx0kyvwp02zppvZrlUrePNN6N7d70hMQ6JV9RwCvAWcGGaaAlbwZ7ncXNdd3rZtVugbx+uulupqd1Jw9dVgGR1ST7Sqnpu9t39S1W+Cp4lI74RGZVJabS08+SSccYbrwNuYUAsXwt//DnvuCRdc4Hc0JlQsT+4+F2bcs2HGmSzx4ouuCd+rr/odiUlJpaUMOrGYL7f25oI/FadNH77ZJFodf39gH6B9SD1/OyC9upQ3cfWzn8Ebb8ARR/gdiUk5Xh++VFVRDLAMFl10B59/3JNT77InvFJFtDr+fsAJQAd2ruffCFycwJhMilq61N3A2203OPJIv6MxKWncuLq0DgE3bhnHh/ftyfG3WfebqSKW7JwHquoHSYonLHty13+qMGoUbNjgcvHkWBc+JpwwffiuowMbaE+RlvkTUxZrTnbOX4tIh6AFdRSRh+IZnEl9IvDgg3DffVbomyjC9LjTkR8pKnLvn3wy5IKgtBSKi+t69bL7AckRy7/wvqr6Y2BAVdcB+yUsIpNyAmkYBgxwT+caE1GUPnw//xzOOgvuv98bH7gfsGyZu0pYtswNhyv87QARV7EU/Dki0jEwICKdiC2rp8kA778P/fq5MzVjGhSlD99Bg+Ctt+CKK7x5w9wPoKrKjQ/WmAOEiUksdfznAjcAz3ijfgGMV9VHExxbHavj98/WrfDXv7oHcSzrpomXdetgQaeDGMV79SeK7NyrV3GxK+xDFRVZt48NiFTHH2ta5gHsSMr2lqouiHN8UVnBn3xbtri/+dZw1yTAGWfAtGd+oKy2F22o3HliaIEe5oYxUP8AYeppzs1dgE5ApareC6yxJ3cz329+Awcf7M74jYm3O+6AKTfMpk1BSIHu3Q/YSZgbxlHHmwY1WPCLyM3AtcD13qg84LFEBmX8N2YMnHKKa7dvTLz17AkH//komDSJObscTzUtdrofsJMoN4xN08Ryxv9z4CRw12Oq+h3QYG2viDwkIhUiMj9o3C0i8q2IzPVexzU1cJMY27e7vyedBNdd528sJvMtGlbCyO9fYeL/bXfVO6GFPkS9YWyaJpaCf5u6GwGBztYLY1z2ZOCYMOMnquoQ7/VajMsySbB6NeyzD7z0kt+RmGzRr597NmTs2AZmLClxB4ba2sgHCBOzWAr+p0XkAaCDiFwMTAP+1dCHVPUd4IdmxmeSrHdv14jCmGQZOxbat3dlemVlw/Ob5muwPb6qThCRI4ENuPw9N6nqm81Y52VeE9HZwNXeA2HGZ6quJ63XX/c7EpONamvhhBOgsBCeftr6d0i0mFr1qOqbqnqNqv6+mYX+/UAfYAiu8/aI/fqKyFgRmS0is9esWdOMVZqGPPYYnH02bN7sdyQmW+XkwFFHwTHhKodN3EVLy/w/VT1IRDbi1e+HWAvcoar/iHVlqloetPx/Aa9EmXcSMAlcO/5Y12Eab/Vq98rL8zsSk83qnug1CRfxjF9VD/L+tlXVdqEvYBjwu8asTER6BA3+HJgfaV6TPL//vcuv38IScZgU8PrrcOyxrktPkxgxVfWIyGARucx77QugqmuB0VE+8wTwAdBPRFaKyIXA/4nI5yIyDzgUuLLZW2CapKYGLroIPvnEDefm+huPMQFbtsB330FFhd+RZK4Gz/FE5He4jlcCnauXisgkVb1HVVdF+pyqnhlm9L+bFqaJt5Ur3Vn+yJEwdKjf0Rizw5gxcPzxdgWaSLF8tRcCB6hqJYCI3I47k78nkYGZxCoqgi++sMRrJjW1aOGqesaPd+lDunXzO6LMEktVjwA1QcM13jiThr78EiZMcM03rdA3qWzJErj9dnjhBb8jyTyxnPH/B/hIRKZ4wz/DqmzS1uTJ8NBDcO65dhZlUtuAAbBoEXW9d5n4iZqWWURygJHAFuAgb/S7qvppEmKrY2mZ40cVVqywxIYmvSxZ4qp+BgzwO5L0Eiktc9QzflWtFZH7VHU/4JOERWcS7vHH4fDDoXt3K/RNeqmpgeOOc7/dd9/1O5rMEEsd/3QROUXEHqJOVxUVLh+KZbE16Sg3Fx5+2J28mPiIpY7/V8BVQI2IeP0yod5DXCntm3XfsGnbJroVdqNLQRdyc7KzsXq3bvDBB7DHHn5HYkzTHHjgjvcbN1rDhOaKJUlb2n7FE96fwD9mu4wSgtCloAvd23SnW2E3uhfu/LdbYbedprXOa+1z9M1XVQUffgiHHQaDBvkdjTHNd+WVMH06zJ4NLVv6HU36iukRCRE5GXdzV3E3d19IZFDxcunwSxldPJryynIqKiso31RORVUFFZUVfPztx1RUVrBx28awn23Tss2Og0Ob7nQrqH9wCAx3yO9AjsTai2Xy3HYb3HorfPWVne2bzHDYYdCxo99RpL8GO1sXkX8AewJPeKNOB75W1d8kOLY6iWzVU7W9ijWVa+oODnUHiMqKHQcM7+/3Vd9Tq/U7d26R04KuBV0jXk0ExgdeLXOTc6pSVQXTprnetIwx2SdSq55YCv4vgb29XrgCTTy/UNW9ExJpGKnSnLOmtoa1m9fudHAId4Ao31ROeWU5W6q3hF1Oh/wOMV9NtG3ZlsbeV//yS3eGb5fCJlN99BH861+uB8ac1LvYThlNas7pWQL0ApZ5wz29cVknNye37qx9YLeBUedVVSq3V9a/eggcMKrc+y8qvuCtyrf4YXP4zspa5baiu7Sh26qNdFu3je4U0m3EoXQfcVi9exNdCrqwubIFo0fD0Ue7lhDGZKKFC12uqRUr7AGvpojljP9tYDjwMa6OfwSu96z1AKqa8IqEVDnjT6TtNdv5vur7egeI8lkzqHh3KuX5NVQUQkUhlLeB7WEaKAlC54LO5C86h133WEfvvTaHrW4KjCtsGWv3ycakFlVXlVkY+hMuLYVx42D5cvfAyvjxWd0/b3PO+G9KQDwmRF5uHj3a9qBH2x47TzjrHlhWs9MoBdb37UnFB9N2uqJYuqKKypZfU1G0jIrKCj5Z5aat37o+7DoL8wp3PigU1K9uCkzv1LpTSt7ANtlJxBX6tbXw6quu20Z5vNQ9sFJV5WZatmxHL+5ZXPiH0+AZfyrIhjP+iHJy3OlNKBH3q/c8/DD89rfw3nv1m25uqd7Cmso1UaucAtPWVK6hRmsIlSu5dC3sWv+gEOEGdn6L/Hh/E8bU89hjcM458OabcMRFxa6wD1VUBGVlyQ4tJTTnjN/4qVev8D/mkLwLhx8OF18cPpdJfot8erbvSc/2PRtcXa3Wsm7zuvpVTpU738z+et3XlG8qp3J7ZdjltGvVLuoBInhah/wOjb6BbQzAGWe4M//DD8dV74QTaXwWszP+VFcacvkKUFDgmjOUlLB5M+TnuwsAP1Ruq2RN1ZrwN7Grdm4e+33V92iY7pvzcvIi3ocIrXLqWtCVvFzrHNjUV9WrPwUrFtWfYGf8TTvjF5HWQC9VDfOtmoQK1E2GuWFVXe16Ktp7b7jvPn/CK2xZSGHLQoo7FDc4b3VtNWur1jb4zMSCNQso31TO1pqtYZfTqXWnBqucAsNtWraxq4ksMG8eHLXhUx5pdTpHbX15x4SCAktSFUYsXS+eCEwAWgK9RWQI8KdktOYxnpKSsDencnLg4IOhTx8fYmqCFjkt6N6mO93bdG9wXlVl47aNUaubKiormFc+j4rKCtZtWRd2Ofkt8mN6ZiLb8zmlu732gtHHtGaXwb+CB+ZZq54GxNKccw5wGDDTS8+MiHyuqknL/pLVVT0R1NbagyvBttVsi3gDO9wBo7q2ut4yIuVzinQTuyCvwIctNSZ2zanq2a6q60Mul1P/xkAG+/xz15LhiSdcNY+Blrkt2a3dbuzWbrcG51VV1m1Z1+DBIZZ8TrEk/OtW2I2OrTtac9gk2LwZbrkFTjkFRozwO5rUFUvB/4WInAXkikhf4HLg/cSGZaKpqnI5yjt08DuS9CQidGrdiU6tO9G/S/8G59+8fXP99BwhB4yl65bywcoPGsznFDY1R8jVRNeCrrRq0SoRm57xtm93efvbtbOCP5pYqnoKgHHAUbhO1l8H/qyq4RPRJIBV9dSn6l9LHhNZcD6naAn/AtM2V28Ou5wO+R1ivppo16qd3cAOsn49tG/vdxSpoclJ2lKBFfzOX/7iOqC4/HIr9DNBaD6ncM1hg6et3bw27HJa5baKeFAIPWB0KehCi5zseHynrMylcM7mg0CT6/hFZC/g90Bx8Pyqelg8AzTR1dbCJ5+4S1iTGUSENi3b0KZTG/p0arhpVmg+p7BXFJvK61o6bavZVn+dXj6nWG5edy/snrb5nCoqYJ994LLL4Pbb/Y4m9cRy6H8G+CfwIFD/WX6TFDk58OyzUF1tZ/vZKmI+pzBUlfVb1zeYQvzT1Z9Svqk8Yj6ngryCsAeIcFcXnVp3SpnmsN26wd13uyy1pr6YmnOq6v5JiiesbK/q+d//XH79XXf1OxKTqbZWb23wBnbgb6R8TjmSU3cDO/i5iXD9TFg+p+RodFWPiHTy3r4sIpcCU4C6RylVNXwCeRNXtbVw7rnQu7fra9SYRGjVolWj8jn9sPmHiFcTgfcfrPiANVVr2LRtU9jltGvVLmKVU+gVRVPzOa1bB5dc4v6Hjjuu0R/PWNGqeubg2usHvu1rgqYpYL24JkFOjus+cVP4/x1jki5HcuhS0IUuBV0Y0DVMVsAQofmcwh0kFv+wmP8t/1/UfE5dC7vWewo73MN23Qq71eVzKiyERYssT1uoWKp68kObboYbl0jZXtVjTLaoqa3h+6rvY3oCO1o+p475HXfcqG69C93bhkkp7k1vSvem6aI5T+6+DwyNYZyJs+eeg1degbvustY8Jjvk5uTW5XMaRPSsMKrKpm2bGjxAzP9+HtPLylm3uB90+Abalu+0nPwW+TE/M5Ep+Zyi1fHvAuwGtBaR/dhR5dMOsCQlSbBihcs62KaN35EYk3pEhLat2tK2VVv27LRn1HnLy6FXL+Xciyq59OrFdVcMdVVQ3jMT3238jrmr51JRWcH22u311+nlc4qU6C/0Bnaq5nOKWNUjIr8EzgOGAbPYUfBvAB5W1eeTESBkd1WPJWMzJj6mToWDDortREpV+XHLjztVK0VqElu+qTzmfE7RUognIp9Tk5/cFZFTVPW5uEbTSNlW8KvCkiXQt6/fkRiTeWpr3bMw8azWD+RzCnRfGi7pX+DgsaZqTcz5nLoVdOOioRexd9emZWNsch2/34V+NnrnHRg9Gl5+2XUibYyJj4oKOOkk90Tv2WfHb7mt81pT1KGIog5FDc5bq7WsrVobNndT8LjFa12V1PF7Hd/kgj+S7EjakWYGDYK//tXrR9QYEzddukD37tC6tX8x5EgOXQu70rWwK/uwT4PzJyKfWrQ6/l+o6jMi0ltVv4n7mhsh26p6jDEmHiJV9US7k3C999eqepJo/Hh47z2/ozAms6m6vP1Ll/odiT+iVfWsFZE3cP3svhQ60frcjb/16+Hee2HbNhg1yu9ojMlcFRUwdiz86lfwt7/5HU3yRSv4j8c9pPUo0OivRkQeAk4AKlR1oDeuE/AULsVzGXCaqobvJTsLtW/vzkBqLAeqMQnVvTt8+CEMaDjjREaKWNWjqttU9UPgJ6r6Ni53zxxVfdsbbshk4JiQcdcB01W1LzDdGza4vkJV3U0ne2DLmMQbONA9I7OtfrcFGS+WpwW6i8inwBfAAhGZIyIDG/qQqr4DhGbwHAM87L1/GPhZI2LNaFddBQcfbGf7xiTTggXueZlp0/yOJLliac45CbhKVWcAiMhob9xPmrC+7qq6ynu/GugeaUYRGQuMBejVq1cTVpVeRo50+fZz0z8NiDFpo08f2H//7LvKjuXJ3c9UdXBD4yJ8thh4JaiO/0dV7RA0fZ2qdmxoOdac0xhjGq8pzTkDlorI/xORYu91I9DURlDlItLDC6gHUNHE5WSMtWvhySetiseYpCotheJiV8lfXEzVQ09y332ua9NsEEvBfwHQFXge16a/izeuKV4Cfum9/yXwYhOXkzEefRTOPNN1FmGMSYLSUteWc9ky16Ji2TKmXfIcl10Gb77pd3DJ0WBVT5MXLPIEMBp3oCgHbgZeAJ4GegHLcM05G+zCMZOremprXbOynzTljokxpvGKi12hH0SBObucwLBVL/sSUqI0OTtnKsjkgt8Yk2Q5Oe5MP5QI1NZSXQ0tMiSLWXPq+E0CbNniWvK8nFknGMakvkitBHv14vHHXfPO9euTG1KyNVjwi0i95AHhxpnGWb3anXhkWzMyY3w3fjwUhPSMVVAA48fTr59r3llZ6U9oyRJLc85PVHVoQ+MSKVOregJffYb282xM6iothXHjYPlydwUwfjyUlPgdVdw1uiMWETkQ95BWVxG5KmhSO8AeM2qGL76AoiI72zfGNyUlUQv6775zWXJ/8YskxpRE0ap6WgJtcAeHtkGvDcCpiQ8tM9XWwqmnul6AjDGp6c9/hvPPz9y6/liqeopUdVnUmRIs06p6PvzQJYb66U/9jsQYE87q1a6ev08fvyNpnib3uQu0EpFJuFTKdfOr6mHxCy+7jBzpdwTGmGh22WXH+5qazMuhFUvB/wzwT+BBwBILNMNbb8H06XD99Va/b0w6uOIKWLECnsuwfghjacdfrar3q+rHqjon8Ep4ZBnovffgscegZUu/IzHGxGL3ik/oPW0SNdLCPfFbWup3SHERrbP1Tt7by3HJ1KYAWwPTY0m1EC+ZVMe/aZOd7RuTFgI5faqqdowrKIBJk9Km6WejUzaIyDe4FBbhWpmrqu4R3xAjy4SCf9066NhgAmpjTMoIyukzl8HkUsMg5ru22GVlvoYWq0bf3FXV3okNKXt8+ikceCBMmQLHHut3NMaYmCxfDsA28jiGqRzIB0zh5Lrx6azBm7sicnKY0euBz1U16/Ppx6JLF7j4YsvAaUxa6dULli2jJdt5npMZwIId49NcLK16LgQOBGZ4w6NxHa/3FpE/qeqjCYotY/TsCffc43cUxphGGT++ro7/J3zgxnk5fdJdLK16WgB7q+opqnoKMABX938AcG0ig8sEDz0EX33ldxTGmEYrKXE3couKQISlux7EEb2X8OmA9LixG00sBX9PVS0PGq7wxv0AbE9MWJlhwwa48kq47z6/IzHGNElJibuRW1tLpy/e5duaHnz7rd9BNV8sVT0zReQV3INcAKd44wqBHxMVWCZo1w4WL7bsm8Zkgg4dYMGCzPh/jqXg/w2usA/k4H8EeE5dO9BDExVYulN1P5Bu3fyOxBgTLyLuf3vuXNhvP7+jaboGq3rUeVZVr/Rez2o69NfosxtugJNPdnk+jDGZY+JEGDYsve/dRcvH/z9VPUhENuJu5tZNwh0P2iU8ujTWubNL9JRpyZ2MyXYlJa7ap3caP+lkna0bY0yGalZn6yJykIic773vIiJpfKxLrMpKePddv6MwxiTa00/Drbf6HUXTxNLZ+s249vrXe6NaAo8lMqh09tBDroOVzz7zOxJjTCLNmAHPPw/V1biEbsXFkJOTFlk8Y2nV83NgP+ATAFX9TkTaJjSqNHbhha5uf/BgvyMxxiTShAnQujXkPBGSxXPZMjcMKZvFM5aqnm1eKx4F8NrvmwgKCjK3g2ZjzA6Fhe4Ef+sNf2RjVUhRWlUF48b5E1gMYin4nxaRB4AOInIxMA34V2LDSj/V1a7AnznT70iMMcmyeTP0X/46f+Km+hNTOItng1U9qjpBRI4ENgD9gJtU9c2ER5Zmli2DTz6BH3/0OxJjTLK0bg1jOzzDiB9frz8xhbN4xlLHj1fQW2EfRZ8+sGiRu/QzxmSP6+/dDcZ+CEEddaV6Fs9oD3CFPrhVNwl7gGsn330H3btDi5gOo8aYjFJSwsbNLXjgmiVc+OPf6FjUzhX6KXpjF6L3wFXXckdEPlXVNM5MkTiq8POfu24Vp071OxpjjB+WDj+da36Ebg+P49xz/Y6mYbGeo6b+470++sMfLDWDMdls8GBX1bvXXn5HEhurnGgmETjlFL+jMMb4LVDoBzLzprJodfzBfe12CO17V1WfT1hUaWL2bPj4Y/fQVqtWfkdjjPHbAw+416xZqV0LEO2M/8Sg92+HDCuQ9QX/M8/Agw/CuedawW+McU/t9+sH69dDp05+RxOZZedsBlVYudJ1pm6MMammWdk5TX3bt7t6PCv0jTGhVq6EJUv8jiIyK/ibYMkSV+BPn+53JMaYVFNd7XrouvZavyOJzFr1NEF1NRx4IOyzj9+RGGNSTYsW8O9/w4ABfkcSWUwFv4j8BCgOnl9VH2nqSkWkDNgI1ADV4eqgUln//jBlit9RGGNS1fHH+x1BdLF0xPIoMAE4CBjuveJRUB+qqkPSrdB/7TVYs8bvKIwxqW7pUrjoIli71u9I6ovljH8YMEDToflPgm3aBKefDqed5i7ljDEmks2b4amnXLr2o4/2O5qdxVLwzwd2AVbFcb0KvCEiCjygqpNCZxCRscBYgF4pkt60TRv3wFbr1n5HYoxJdfvsA6tWuXIj1cRS8HcBFojIx8DWwEhVPakZ6z1IVb8VkW7AmyLypaq+EzyDdzCYBK4dfzPWFVd77+13BMaYdBEo9CsrXY9dqSKWgv+WeK9UVb/1/laIyBRgBPBO9E/5a+JE14zzrrss/bIxJnZXXglvvAGff546/XXE0gPX2/Fcoddnb46qbvTeHwX8KZ7rSISKCvdQhhX6xpjGOOQQ6NLFPfSZKqldGkzZICIjgXuAvYGWQC5Q2dSOWERkDyDQGLIF8LiqRu2qJlVSNqRD1j1jjAmIlLIhlvPXe4EzgGdwLXzOBZqcdVpVlwKDm/r5ZNu2zfWn27evFfrGmKZRhRkzYNdd3XNAfoupxklVlwC5qlqjqv8BjklsWKnjscfcjpo71+9IjDHpauNGGDMG7vzNYigudpX9xcVQWupLPLGc8VeJSEtgroj8H65ZZ4rcoki844+HCRNcDzvGGNMU7drBm7+fypDbz4LN69zIZctg7Fj3Psn988ZSx18ElOPq968E2gP/8K4CkiJV6viNMabJiotdYR+qqAjKyhKyyianZVbVZYAAPVT1j6p6VTILfb+owlVXwaef+h2JMSYjLF/O/xjFobzFBtruND7ZYsnVcyIwF5jqDQ8RkZcSHJfvli6FyZNh3jy/IzHGZIRevWjFVlbRg2UU7TQ+2WJ9gGsEMBNAVeeKSO8ExpQS+vRxV2X5+X5HYozJCOPHM3zsWBZUDSAHr4q9oADGR23NnhCx3KTdrqrrQ8alTAqFRNi0yf1t2xby8vyNxRiTIUpKYNIkcop6UU0LvtttOEyalPQbuxDbGf8XInIWkCsifYHLgfcTG5a/Tj4Z2rd3nakbY0zclJRASQmH/RRyc2FG8st8ILaC/7fAOFyCtieA14E/JzIoP6nCiSdaBk5jTOL89rfQsqV/62+wOWcqsOacxhjTeI1O2dBQy51mpmVOSV99BQsXujP+VMmiZ4zJTBs2wCOPwLnnuge8kilaVc+BwApc9c5HuLb8Ge3++929luXLoXNnv6MxxmSyRYtclU/nznDmmcldd8SqHhHJBY4EzgT2BV4FnlDVL5IXnpOsqp7qatduf+jQhK/KGGOYPx8GDkzc8hv95K6XkG2qqv4SGAksAWaKyGWJC9M/qi7XvhX6xphkSWShH03UmmwRaSUiJwOPAb8B7mZHLv2MsWqVS8L23nt+R2KMyTa33QYXX5zcdUa7ufsIMBB4Dfijqs5PWlRJVl7untDt0cPvSIwx2aay0t3ora1NXqOSaHX8tUClNxg8kwDa1B64msKacxpjMlUie/ZrSh1/jqq29V7tgl5tk1noJ9rs2bBli99RGGOyVaDQLy+HrVuTs86sbq1eWQlHHQW//rXfkRhjstn8+dCzZ/LSxMSSsiFjFRTAc89B165+R2KMyWYDBsC4cTByZHLWl9UFvwgceqjfURhjsl1ODtx8cxLXl7xVpZZHH3VpsKur/Y7EGGOczz93ZVOiZW3B//778OqrLjWqMcakgvvug9/9LvENTrI6O+fmzZZ+2RiTOlatcuma45UrrNHZOTNVbS2sW+e+WCv0jTGpJFkPkWZdVc8LL7i+jT/91O9IjDGmvtWrXWr4115L3DqyruAfONC12993X78jMcaY+jp3hm+/dTUTiZJ1VT177QV/+5vfURhjTHh5eTBnTuLSOECWnfH//e9QVuZ3FMYYE52Iy+GzfHlilp81Bf/SpXDddfDii35HYowxDbvuOpcuvqoq/svOmqqePfZwhX+nTn5HYowxDTvtNOjTJzGpmrOi4A/kud59d78jMcaY2Oy/v3slQlZU9ZSUwK9+5XcUxhiTGjL+jF/VVfMUFvodiTHGpIaML/hFXDI2Y4wxTkZX9Xz7Lcya5XcUxhiTWjK64F+xAi68ECoq/I7EGGNSR0ZX9YwcCZ99ltgn4IwxJt1k9Bk/WKFvjDGhfCn4ReQYEVkkIktE5LqErKS0FIqLXQP+4mI3bIwx6SKBZVjSq3pEJBe4DzgSWAnMEpGXVHVB3FZSWgpjx+541nnZMjcMrlG/McaksgSXYX6c8Y8AlqjqUlXdBjwJjInrGsaNq5/goqrKjTfGmFSX4DLMj4J/N2BF0PBKb9xORGSsiMwWkdlr1qxp3BoipbRLVKo7Y4yJpwSXYSl7c1dVJ6nqMFUd1rVr18Z9uFevxo03xphUkuAyzI+C/1ugZ9Dw7t64+Bk/HgoKdh5XUGCP8Bpj0kOCyzA/Cv5ZQF8R6S0iLYEzgJfiuoaSEpg0CYqKXHvOoiI3bDd2jTHpIMFlmKhqXBbUqJWKHAfcCeQCD6lq1MPYsGHDdPbs2ckIzRhjMoaIzFHVYaHjfXlyV1VfAxLYh7wxxphIUvbmrjHGmMSwgt8YY7KMFfzGGJNlrOA3xpgs40urnsYSkTXAsiZ+vAvwfRzDSUWZvo22fekv07cxVbevSFXrPQGbFgV/c4jI7HDNmTJJpm+jbV/6y/RtTLfts6oeY4zJMlbwG2NMlsmGgn+S3wEkQaZvo21f+sv0bUyr7cv4On5jjDE7y4YzfmOMMUGs4DfGmCyT0QV/Ujp1TzAR6SkiM0RkgYh8ISK/88Z3EpE3RWSx97ejN15E5G5vm+eJyFB/tyA2IpIrIp+KyCvecG8R+cjbjqe8FN6ISCtveIk3vdjXwGMkIh1E5FkR+VJEForIgZm0D0XkSu/3OV9EnhCR/HTfhyLykIhUiMj8oHGN3mci8ktv/sUi8ks/tiVUxhb8QZ26HwsMAM4UkQH+RtUk1cDVqjoAGAn8xtuO64DpqtoXmO4Ng9vevt5rLHB/8kNukt8BC4OGbwcmquqewDrgQm/8hcA6b/xEb750cBcwVVX7A4Nx25oR+1BEdgMuB4ap6kBcuvUzSP99OBk4JmRco/aZiHQCbgYOwPU3fnPgYOErVc3IF3Ag8HrQ8PXA9X7HFYftehE4ElgE9PDG9QAWee8fAM4Mmr9uvlR94Xphmw4cBrwCCO4pyBah+xJ4HTjQe9/Cm0/83oYGtq898E1onJmyD9nRj3Ynb5+8AhydCfsQKAbmN3WfAWcCDwSN32k+v14Ze8ZPjJ26pxPvkng/4COgu6qu8iatBrp779Nxu+8E/gDUesOdgR9VtdobDt6Guu3zpq/35k9lvYE1wH+86qwHRaSQDNmHqvotMAFYDqzC7ZM5ZNY+DGjsPkvJfZnJBX9GEZE2wHPAFaq6IXiaulOJtGyXKyInABWqOsfvWBKoBTAUuF9V9wMq2VFFAKT9PuwIjMEd4HYFCqlfRZJx0nmfZXLBn/hO3ZNERPJwhX6pqj7vjS4XkR7e9B5AhTc+3bZ7FHCSiJQBT+Kqe+4COohIoIe44G2o2z5ventgbTIDboKVwEpV/cgbfhZ3IMiUfXgE8I2qrlHV7cDzuP2aSfswoLH7LCX3ZSYX/Inv1D0JRESAfwMLVfXvQZNeAgItBH6Jq/sPjD/Xa2UwElgfdGmaclT1elXdXVWLcfvoLVUtAWYAp3qzhW5fYLtP9eZP6bMuVV0NrBCRft6ow4EFZMg+xFXxjBSRAu/3Gti+jNmHQRq7z14HjhKRjt6V0VHeOH/5fZMhkS/gOOAr4GtgnN/xNHEbDsJdTs4D5nqv43B1otOBxcA0oJM3v+BaM30NfI5raeH7dsS4raOBV7z3ewAfA0uAZ4BW3vh8b3iJN30Pv+OOcduGALO9/fgC0DGT9iHwR+BLYD7wKNAq3fch8ATunsV23FXbhU3ZZ8AF3rYuAc73e7tU1VI2GGNMtsnkqh5jjDFhWMFvjDFZxgp+Y4zJMlbwG2NMlrGC3xhjsowV/CYuRGRTEz93nojcG+dYWonINBGZKyKnx3PZIespDs7cmCyxfmciMlNcdtrPROS9oOcIws27q4g8G8Myb2hsvCb1WMFvMtF+AKo6RFWfCp7gZW3NJiWqOhh4GLgj0kyq+p2qnhppehAr+DOAFfwmrkRktHemGcg9X+o9zYmIDBeR970z0I9FpK33sV1FZKqXr/z/gpZ1lIh8ICKfiMgzXr4iROQ2cf0TzBORCSHr7wY8Bgz3zvj7iEiZiNwuIp8AvxCRM0Xkc3G5428P+uwmEblDXF75aSIywtuWpSJyUiO+g4tFZJa3nc+JSIE3frKInBo036Z4f2dRvAPs6T1Zeoe37Z8HroiCr168K4rnQ5cvIrcBrb3vtTTW78OkIL+fILNXZryATd7f0bhsi7vjTiw+wD193BJYCgz35muHS152nje+Pe6JzmW43CZdcIVVoTf/tcBNuCcnF7Gjv+gOYWIZjfcEsDdcBvzBe78rLsVAV2/9bwE/86YpcKz3fgrwBpCHy58/N8x6iglK2Rs0vnPQ+78Av/XeTwZOTdR3FiaOmXhPkALXAE8BpwBv4nLmd/e+ix7B2xJt+YGY7ZXer0ACJWPi6WNVXQkgInNxhcp6YJWqzgJQL8Ood2I7XVXXe8MLgCKgA64Dnfe8eVriCsT1wBbg3+J663olxpgCVT7DgZmqusZbXynwU1wahW3AVG++z4GtqrpdRD73tiFWA0XkL942tCG23Czx+M5W1FsqlIrIZtzB77fAVcATqlqDSzj2Nu47mRfyuViXb9KQFfwmEbYGva+h4d9ZuPkFeFNVzwydWURG4BKBnQpchsvo2ZDKGObZrt5pLa5vgK0AqlorO7JMxmIy7iriMxE5D3dGD643tRwAEcnBHcwC4vGdhVOiqrMDA95BIxaNjcekEavjN8myCOghIsMBRKRtA4Xph8AoEdnTm79QRPby6vnbq+prwJW4apjG+Bg4RES6eDd6zwTebuzGNKAtsEpcOu2SoPFlwP7e+5Nw1UjRNPY7i8W7wOni+jjuirva+bgRn9/ubZdJY3YUN0mhqtu8G4n3iEhrYDMuj3uk+dd4Z8tPiEgrb/SNwEbgRRHJx10VXNXIOFaJyHW4lMECvKqqLzbwsWj6icjKoOErgf+H6yVtjfc3cEP2X17sn+GqlKJehTT2O4vRFFw3iJ/h7mn8QVVXS+wdnk8C5onIJ+rSZ5s0ZNk5jTEmy1hVjzHGZBkr+I0xJstYwW+MMVnGCn5jjMkyVvAbY0yWsYLfGGOyjBX8xhiTZf4/leGnWB7BwnMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "processTrajectories(d2hs,trials)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A quick look at the plot on the right makes it quite clear that\n",
    "\n",
    "* a **quadratic** fit is far better than a linear one\n",
    "\n",
    "<b style=\"color:blue;font-size:120%\">how bad a fit is the line and how good is the quadratic fit?</b>\n",
    "\n",
    "* **Coefficient of Determination(决定系数)**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "####  something about plotting lines with matplotlib.pyplot\n",
    "\n",
    "* [Unit2-2-PLOTTING-USING-MATPLOTLIB](./Unit2-2-PLOTTING-USING-MATPLOTLIB.ipynb)\n",
    "\n",
    "The function `processTrajectories` plot lines with a legend(图例说明)  \n",
    "\n",
    "**line styles**\n",
    "\n",
    "* the linear fit for the data.:  plt.plot(distances, altitudes, 'g', label = 'Linear Fit')\n",
    "\n",
    "  *  <b style=\"color:blue\">green solid line</b> \n",
    "\n",
    "* the quadratic fit for the data.:  plt.plot(distances, altitudes, 'b:', label = 'Quadratic Fit')\n",
    "\n",
    "   * <b style=\"color:blue\">blue dot line</b> \n",
    "   \n",
    "**Legend of matplotlib.pyplot**\n",
    "\n",
    "* **label**：a legend of line      \n",
    "       \n",
    "* **Place a legend on the axes**\n",
    "\n",
    "  *  `plt.legend(loc ='best')` # Place a legend on the axes.\n",
    "\n",
    "Location String： **`'best'`**   \n",
    "\n",
    "> reference: [matplotlib.pyplot.legend](https://matplotlib.org/devdocs/api/_as_gen/matplotlib.pyplot.legend.html): \n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3 Coefficient of Determination\n",
    "\n",
    "\n",
    "\n",
    "Asking about **the goodness of a fit** is equivalent to asking about the accuracy of these predictions.\n",
    "\n",
    "<b>The mean square error</b> is useful for comparing the <b>relative goodness</b> of two fits to the same data, \n",
    "\n",
    "it is <b>not particularly useful for getting a sense of the absolute goodness of a fit</b>.\n",
    "\n",
    "We can calculate <b>the absolute goodness of a fit</b> using the <b>coefficient of determination(确定系数)</b>, often written as\n",
    "\n",
    "$R^2$\n",
    "\n",
    "Let:\n",
    "\n",
    "* $y_i$ be the $i^{th}$ observed value,\n",
    "\n",
    "* $p_i$ be the corresponding value predicted by model, and \n",
    "\n",
    "* $\\mu$ be the **mean** of the observed values.\n",
    "\n",
    "$$R^2=1-\\frac{\\sum_{i}(y_i-p_i)^2}{\\sum_{i}(y_i-\\mu)^2}$$\n",
    "\n",
    "By comparing\n",
    "\n",
    "<b>the estimation errors</b> (<b>RSS: residual sum of squares</b>) \n",
    "\n",
    "$\\sum_{i}(y_i-p_i)^2$\n",
    "\n",
    "with\n",
    "\n",
    "<b>the variability of the original values</b> ( <b>TSS: total sum of squares</b>), \n",
    "\n",
    "$\\sum_{i}(y_i-\\mu)^2$\n",
    "\n",
    "$R^2$ is intended to capture <b>the proportion of variability </b> in a data set that is accounted for by the statistical model provided by the fit.\n",
    "\n",
    "Its <b>compactness</b> stems from the expressiveness of the operations on <b>arrays</b>."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def rSquared(measured, predicted):\n",
    "    \"\"\"Assumes measured a one-dimensional array of measured values\n",
    "               predicted a one-dimensional array of predicted values\n",
    "       Returns coefficient of determination\"\"\"\n",
    "    # RSS: residual sum of squares\n",
    "    estimateError = ((predicted - measured)**2).sum()\n",
    "    \n",
    "    # TSS: total sum of squares\n",
    "    meanOfMeasured = sum(measured)/float(len(measured))\n",
    "    variability = ((measured - meanOfMeasured)**2).sum()\n",
    "    \n",
    "    return 1 - estimateError/variability"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def processTrajectories(d2hs,trials):\n",
    "    numTrials = trials\n",
    "    distances=[]\n",
    "    meanHeights=[]\n",
    "    \n",
    "    for item in d2hs:\n",
    "        distances.append(item['d'])\n",
    "        meanHeights.append(np.mean(item['h']))\n",
    "    \n",
    "    plt.title('Trajectory of Projectile (Mean Height of {} Trials'.format(trials))\n",
    "    plt.xlabel('Inches from Launch Point')\n",
    "    plt.ylabel('Mean Height of the projectile')\n",
    "   \n",
    "    # the experimental data\n",
    "    plt.plot( distances, meanHeights, 'ro')\n",
    "   \n",
    "    # Linear Fit\n",
    "    a,b = np.polyfit(distances, meanHeights, 1)\n",
    "    altitudes = a*np.array(distances) + b\n",
    "    \n",
    "    r1line=rSquared(meanHeights, altitudes)\n",
    "    plt.plot(distances, altitudes, 'g', label = \"Linear Fit, $R^2$= {:.2f}\".format(r1line))\n",
    "    \n",
    "    # Linear Fit\n",
    "    a,b,c = np.polyfit(distances, meanHeights, 2)\n",
    " \n",
    "    altitudes = a*(np.array(distances)**2) +  b*np.array(distances) + c\n",
    "             \n",
    "    r2quad=rSquared(meanHeights, altitudes)          \n",
    "    plt.plot(distances, altitudes, 'b:',label =\"Quadratic Fit, $R^2$= {:.2f}\".format(r2quad))\n",
    "    plt.legend()\n",
    "    plt.show()\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEWCAYAAABhffzLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAABFFElEQVR4nO2deZgU1bXAf2fYZ1iGHURgEBFZFERB3FGUKEYx4j6uQTGJxojGuJCnPhMSjUTFNaIhiuC+L3kaAVdECYgKgsi+wwCyzgzLMOf9cauHnp7unp6lu3o5v++rr+veulX31NKnbp177zmiqhiGYRiZQ5bfAhiGYRiJxRS/YRhGhmGK3zAMI8MwxW8YhpFhmOI3DMPIMEzxG4ZhZBim+OOIiPyfiFzhtxy1iYj8WkQ2iMhOEWkZx3pOEJGF8T6uiCwXkVOreawGIjJfRNrXnoTJhYj8Q0T+J8ayz4jIn+Moyy9EZJX37B0Rr3rC1LtTRA6KoVyeiKiI1E2EXDXBFH8I3k0OLKUiUhyUzq/KsVT1DFV9toby3C0ik2pyjNpCROoBDwBDVLWxqm4O2R548APXa7mI3FadulT1M1XtXgsyq4gcXNvH9RgJfKqq67y6nvHqGxYiw4Ne/pW1VG/MhJ6/lxfzM6Wqv1LVP8VLlioyFrjee/bmRKnnJK+usC8hr0EWeEb3isieoPQ/Qst79S2tgdxJR9K/mRKNqjYOrIvIcuBqVZ0SWk5E6qpqSSJlqw61LGdboCHwfSXlclW1RESOAaaKyDeq+n4c5fKLX+GUfzA/ApcDb4E7T+ACYEliRUtLOlPJs+c1TsYBX0Uqo6pnBJV/Blitqn8Mc6x0eEbDYi3+GBGRQSKyWkRuFZH1wL9EpLmIvCsiG0Vki7d+YNA+H4vI1UHpX4rIAq/sByLSOWhbLxH5UER+8kwpd4jI6cAdwIVea+Rbr+wBIvK2V3axiFwTdJy7ReRVEZkkItuB20SkKNgsIyL9PJnrhTnPBiLykIis9ZaHvLxDgICJZKuITKvsmqnqDNwftXeE6xe2ruDrHSTXASLymif3MhG5IWhbHe96LRGRHSIyW0Q6isinXpFvvet3YehxQ849S0Ru846zWUReFpEWEcp2Ag6iooJ5BzheRJp76dOB74D1IftHexbGiTNpbPfO5YSgbXd7ck30zvV7ETkqym2oFBE5NOjZWygiFwRtK2e+EZE/iMg6735dLRVb8c1F5D1Ptq9EpKu3X4V7EUaOLBH5o4isEJEC7xybec/JTqCOt3+0l+jNwH+AH6p5LVRErhORRcCioLyDvfUzRWSOd29WicjdUY51pYgs9a7FMqmixSCemOKvGu2AFriWx0jc9fuXl+4EFAOPhttR3Of/HcC5QGvgM+AFb1sTYArwPnAAcDAw1Wsl/wV4yfvc7OMd7kVgtVf2POAvInJKUHXDgFeBXODvwMe4VmeAy4AXVXVvGFFHAwOBvkAfYADwR1X9EejllclV1VPC7Bt8viIix3n7BD7LQ69f2LrCHCsLp1C/BToAg4EbReRnXpGbgIuBoUBT4JdAkaqe6G3v412/l6LJDPwWOAc4CXdttwCPRSh7GLA0TItwF661f5GXvhyYGHI+EZ8Fj//irkkL4HngFRFpGLT9bNwzkAu8TYRnLhZEJAf40KunjSf34yLSM0zZ03HX+lTcMzoozCEvAv4XaA4sBsYAxHgvrvSWk3Ev1cbAo6q6O+hLvI+qdo1wLp1x9/6eqCddOecARwMVrgFQiLunucCZwK9F5JwwsuQADwNnqGoT4FjgmxrKVXuoqi0RFmA5cKq3PgjYAzSMUr4vsCUo/THOVATwf8CIoG1ZQBFOCV4MzIlwzLuBSUHpjsA+oElQ3l+BZ4LKfxpyjAuB6d56HVzrc0CE+pYAQ4PSPwOWe+t5gAJ1I+wb2L4VpzQXADdEun6V1DUI9wkO7k+4MqSu24F/eesLgWERZFLg4KB02XHD3OMFwOCgbe2BveHOF8gHvgzJewb4M3A8MAOnHDYAjYDPgSsrexYinMMWnMIL3N8pQdt6AsVRnkkFtnv3JLDsCjxT3rPxWcg+TwJ3BZ+Ttz4B+GtQuYODr69X9umg7UOBHyLdizCyTgV+E5TuHnz9Y9j/LeDCULmjLaHlvDpOifYMhWx7CHgw9P8B5HjXejjQqDI5Er1Yi79qbFTVXYGEiGSLyJPep+l24FMgV0TqhNm3MzBORLaKyFbgJ0BwLdiOxG4DPgD4SVV3BOWt8I4TYFXIPm8BPUWkC3AasE1VZ0Y5/oqQYx8Qo2wBWqlqc1XtoaoPB+WXu35VqKszcEDg2nnX7w5cnwNU7fpFozPwRlAdC3Av2bZhym4BmoQ7iKp+jmvJjwbeVdXiMPVEehYQkd97ZqBt3vZmQKug/YPNRkVAQ4k+kqSfquYGFuDeEFmODrm2+bivs1AOoPyzFfqchZOtcZgykQj3PNQl/PUvh4ichWsMVfZVFwvhzitQz9Ei8pE4k+M2XD9Pq9ByqlqIe6n+Cljnmb8OrQXZagXr3K0aoa5Mb8a1So5W1fUi0hdn1pAw+64Cxqjq5NAN3ifqRRV3CVvnWqCFiDQJUv6dgDWR9lHVXSLyMnApcCjwXIS6AscP7kTr5OXVBuHOJZa6VgHLVLVbhOOuAroC82oo3yrgl6o6PYay3wFdJHIH4CTgTpzZIlw9kZ6FE4A/4MxZ36tqqYhsIfwzVRusAj5R1dNiKLsOODAo3bGWZQk8DwE6ASW4r6bKGAwcJa7/CNzLcp+IHKaqw6LsF45oLoufx5nWzvD+Vw8RRvEDqOoHwAci0gj3JfgUcEK4sonGWvw1ownOrr9VXCfgXVHK/gO4XUR6AXidVud7294F2ovIjV5HVhMROdrbtgHI8+zcqOoq4AvgryLSUEQOB0bgFE00JuLsp2cTXfG/APxRRFqLSCuc8orXcNJY65oJ7BDXMdxIXGdubxHp721/GviTiHTz+hYOl/2d2Rtw9uJY+AcwxnsR48kVVmmo6mqcDXtAhGM9jPu6+jTMtmjPQhOcstsI1BWRO3H9FvHiXeAQEblMROp5S38R6RGm7MvAVSLSQ0SygZjG9wdR2b14ARglIl1EpDH7+7diGVnzP8AhOHNrX1zfx1PAVVWUsTKa4L64d4nIAOCScIVEpK2IDPNs/buBnUBpLctSbUzx14yHcPbbTcCXuM7ZsKjqG8B9wIueWWgecIa3bQdOSZyF+1RexP6W4ive72YR+dpbvxhnT1wLvIGzx1YYchpS/3Tcg/e1qq6IUvTPwCxci3Yu8LWXFw9iqktV9wE/x/2hl+Gu99O4Vh24uQUv40ZzbAf+ibsv4Gziz3pmjOAO7nCMwymM/4jIDtw9PTpK+SdxHeUVUNWfVHWqesbfkG0RnwXgA9xz9CPO1LGLKKaHmuI9e0NwX5xrcc/ffUCDMGX/D/dC+wj30vvS27Q7xuruJvq9mIBrlHyKu8+7cB3uMZ2Hqq4PLLgGWaGq/hSjbLHyG+Ae7/m4E/fchSML1xG+FmfKOwn4dS3LUm0kzHNp1BLihrA9raoTKy2cAMQNwXxeVZ/2W5bK8EYpPa2qsbbWE464oadzcB3C6/yWJ9F4XwXzgAYxtsqNJMFa/HHC+xQ+CNdy8R3PLNIPqI3Or0TQmyS5dpFQN8ywZyYpfXFuExqIm6dwH/COKf3UwxR/HBCRNrhP5k9ww/h8RUSexc0TuDFkNFBSIiLjgFG48eBGcnEtUIAbRbWPJDJfGLFjph7DMIwMw1r8hmEYGUZKjONv1aqV5uXl+S2GYRhGSjF79uxNqto6ND8lFH9eXh6zZs3yWwzDMIyUQkTCDt02U49hGEaGYYrfMAwjwzDFbxiGkWGY4jcMw8gwTPEbhmFkGKb4DcMwMgxT/IZhGBmGKX7DqCazZsEjj8CePX5LYhhVwxS/YVSB0qBQGk8/DY8+CvXqufT8+bA3XPh6w0gyTPEbRoxMnw5dusCPP7r0PffA55+DCJSUwKmnwpVX+iqiYcSEKX7DiML338OCBW794IOhRw8oKnLpNm2gtecFRcR9AfzWixf1008wZAjMjBTS3jB8JCV89RiGH+zZAyedBKecAi+/DG3bwvsRgmvWqQNDh+5PL10KS5ZA/fouvWoVbNsGvXvHX27DqAxr8RtGEC++CJde6tbr14dXX4XHH6/6cY46ChYvhr59XXrcODjiCPclYBh+EzfFLyIdReQjEZkvIt+LyO+8/BYi8qGILPJ+m8dLBsOIheXLYd8+t75xIyxa5FrnAIMGQatW1TuuyP71226D116DFi1c+qqr4I47qiuxYdSMeLb4S4CbVbUnMBC4TkR6ArcBU1W1GzDVSxuGL3zxBXTtCu+849K/+Q189RU0a1a79bRqBWef7dZV3UigOnX2b3/xxf0vG8OIN3FT/Kq6TlW/9tZ3AAuADsAw4Fmv2LPAOfGSwTBCKS2FZ5+F11936QED3Oic/v1dOlgZxwsRGD8e/vQnl54/Hy6+GCZO3C+jRUQ14klCbPwikgccAXwFtFXVdd6m9UDbCPuMFJFZIjJr48aNiRDTSGN273a/IvDYY/uVbN26MHo0dOjgn2w9e7rJYJdd5tLvvAO9esGyZf7JZKQ3cVf8ItIYeA24UVW3B29TF+k9bNtGVcer6lGqelTr1hUihxlGzDzyiDPn7NrlFP+778Ibb/gtVXmOPBJyc916Tg506wYdO7r0++/DlCn2FWDUHnFV/CJSD6f0J6uq93HNBhFp721vDxTEUwYj8ygpcR2pgQ/FPn3g3HPLj78P7nhNNk49Fd56y32NAPzlL64jOCBzcbF/shnpQTxH9QjwT2CBqj4QtOlt4Apv/QrgrXjJYGQmS5bAeefBpEkufeKJ8PDD+0fUpBr/+Q+88IJb373bzR6+/35/ZTJSm3i2+I8DLgNOEZFvvGUocC9wmogsAk710oZRI26+GW65xa137w6ffQY33OCvTLVFw4bOVAXOXHXlla5TGmDDBrjzTvdrGLESt5m7qvo5EOmDenC86jUyA1WYMwf69XPpoqL9phGA44/3R65406wZ3BvUVPr4YxgzBi65xM0s3rLF9REEZgwbRjhs5q6Rkvz972527NKlLv34464TN9O48EJYswYOPdSlR4+GQw4xL6FGdEzxGynBli1w663w3/+69CWXuCGZBx7o0gnvrJ08GfLyICvL/U6enGAB9tOu3f71c8+Fm27a7yr6j390HcWGEYwpfiNpUd0/m7VePfjnP91MW4ADDnA+dXwxaUyeDCNHwooVTsgVK1zaR+Uf4NRT9/dt7NrlnMsFXpaq8MMP/slmJA+m+I3EUoWW8gUX7Hdz0Lix86nzu98lQshKGD16/9jQAEVFLj+JaNjQKfqAWF995dxKv/aav3IZ/mOK30gclbSUCwqcF8tAlKthw+Cii/ZPXGrc2Ce5Q1m5smr5PpKVBY0aufVDDoEHHnBxAgDefBOuvRa2b4+4u5GmmOI3EkeElrLe4ZqkU6fCjTfuN01cein8+tdJONmqU6eq5ScJLVrAqFHQpIlLL1nioooFXqj2AsgcTPEbiSOkRbyFXE5mGv9ceRrgJl0tWABHH+2HcFVgzBjIzi6fl53t8lOIm2+Gb791XwWbNrl4AX/7m99SGYnAFL+RODp1QoHZuMH3uWylKdtp2DIHcB24gWGJSU1+vnOv2bmz+xzp3Nml8/P9lqzK1KkDTJ5M7pFdGbr0UQY9OCwpOqmN+CKaAp6fjjrqKJ01a5bfYhg1ZfJk7r9qPn/f+1tW0JkG7HEt5RRVmmlBoN8l2ASXnc27v36PQXcPSp5+FaNaiMhsVT0qNN9a/EbiyM9n2F+OZkTTV6nP3pRuKacNYfpdVhW14Ny/H8s99/gkkxF3rMVvJITi4v2jS4wkIisrrL/njziZATunkZPjg0xGrWEtfsM3iorg2GPh7rv9lsSoQISRSCd3XkpODuzZ46KDzZyZYLmMuGKK34g79eo518gBj5JGElHJCKUNG5zSX7jQB9mMuGGK34grpaVO8Y8bB0OH+i2NUYFKRih17Ajz5nlhISdPpqDjkUnhn8ioGab4jbjx2WeulZ+EE1qNYPLznT+M0lL3G9LZ3qgRMHkyi66+j+6rp/APHZlU/omMqlOp4heRbBH5HxF5ykt3E5Gfx180I9XZs8c5UWvWzG9JjBozejR5uxZwFf9iKP92eUnon8iIjUpH9YjIS8Bs4HJV7S0i2cAXqto3AfIBNqonlVFNQpcLRtUJGf3jJuIdyVHy9X7nSkbSUZNRPV1V9W/AXgBVLSJyZC3DYNQoeOYZt25KP00IGf3zHJfRn1l82uY8nwQyakIsoRf3iEgj3EseEekK7I6rVEbKsmsXfPedhf5LO8aMKTfD9wJeZke9lhx//zCfBTOqQyyK/y7gfaCjiEzGBVG/Mp5CGalLw4bwwQd+S2HUOoEO39GjYeVKGnZqx3VjjoL8fDZtgvnz3ZBdIzWIaeauiLQEBuJMPF+q6qZ4CxaM2fiTn40b4a674L779rv9NTKD/Hx4/31YtgyaNvVbGiOYKtv4RaRfYAE6A+uAtUAnL88wyvjkExcDd8kSvyUxEs1DD7m4vqb0U4dopp6/R9mmwCm1LIuRwpx3Hpx0ErRu7bckRqJp3Xr/ff+//3OTvnr39lcmIzoRFb+qnpxIQYzU5M03oV07GDjQlH6ms3u3i5h22GHwzjt+S2NEI6LiF5FTVHWaiJwbbruqvh4/sYxUYN8+uPNOp/CnTLGhm5lOgwbw4YfQtq3fkhiVEc3UcxIwDTgrzDYFTPFnOHXqwMcfuxm6pvQNgG7d3G9JiWsU3HwztGzpr0xGRaKZeu7yVu9R1WXB20SkS1ylMpKa0lJ48UW46CIXwNswQlmwAB54AA4+GH75S7+lMUKJZebua2HyXq1tQYzU4a233BC+997zWxIjKZk8mcPOyuOH3V345T155sgtCYlm4z8U6AU0C7HzNwUaxlswI3k55xz4z3/g1FP9lsRIOoJi+OYBrICFV9/P3JkdOW+czfBKFqLZ+LsDPwdyKW/n3wFcE0eZjCRl6VLXgdehA5x2mt/SGElJmBi+f9w1mi8fO5gz77Xwm8lCLN45j1HVGQmSJyw2c9d/VOG442D7dueLJ8siORjhCBPDdwu5bKcZnXW5PzJlMDXxzvkrEckNOlBzEZlQm8IZyY8IPP00PPaYKX0jCmFi+DZnK507u/UXXwz5IJg82UXzsqheCSWWv/Dhqro1kFDVLcARcZPISDoCbhh69nSzcw0jIlFi+M6dC5dcAk884eUH+gNWrHBfCdGietkLolaJRfFniUjzQEJEWhCbV08jDfjiC+je3bXUDKNSosTwPewwmDYNbrzRKxumPyBsVK+qvCCMmIjFxn85cAfwipd1PjBGVZ+Ls2xlmI3fP3bvhr/+1U3EMa+bRm2xZQvMb3E8xzG94kaR8lG98vKcsg+lc2cXI9iISCQbf6xumXuy3ynbNFWdX8vyRcUUf+LZtcv9NrSBu0YcuOgimPLKTywv7URjCstvDFXoYTqMgYovCKMCNencBWgBFKrqo8BGm7mb/lx3HZxwgmvxG0Ztc//98MYds2icHaLQvf6AcoTpMI6ab1RKpYpfRO4CbgVu97LqAZPiKZThP8OGwfDhbty+YdQ2HTvCCX8aAuPHM7vdmZRQt1x/QDmidBgb1SOWFv8vgLPBfY+p6lqgUmuviEwQkQIRmReUd7eIrBGRb7xlaHUFN+LD3r3u9+yz4bbb/JXFSH8WHpXPwE3v8uDf9jrzTqjSh6gdxkb1iEXx71HXERAItp4T47GfAU4Pk/+gqvb1ln/HeCwjAaxfD716wdtv+y2JkSl07+7mhowcWUnB/Hz3YigtjfyCMGImFsX/sog8CeSKyDXAFOCpynZS1U+Bn2oon5FgunRxgygMI1GMHAnNmjmdXlhYeXmj5lQ6Hl9Vx4rIacB2nP+eO1X1wxrUeb03RHQWcLM3IczwGVUXSeuDD/yWxMhESkvh5z+HnBx4+WWL7xBvYhrVo6ofquotqvr7Gir9J4CuQF9c8PaIcX1FZKSIzBKRWRs3bqxBlUZlTJoEl14KxcV+S2JkKllZMGQInB7OOGzUOtHcMn+uqseLyA48+34Im4H7VfXxWCtT1Q1Bx38KeDdK2fHAeHDj+GOtw6g669e7pV49vyUxMpmyGb1G3InY4lfV473fJqraNHQBjgJ+V5XKRKR9UPIXwLxIZY3E8fvfO//6dc0Rh5EEfPABnHGGC+lpxIeYTD0i0kdErveWwwFUdTMwKMo+LwAzgO4islpERgB/E5G5IvIdcDIwqsZnYFSLffvg6qvh669duk4df+UxjAC7dsHatVBQ4Lck6UulbTwR+R0u8EoguPpkERmvqo+o6rpI+6nqxWGy/1k9MY3aZvVq18ofOBD69fNbGsPYz7BhcOaZ9gUaT2K5tCOAo1W1EEBE7sO15B+Jp2BGfOncGb7/3hyvGclJ3brO1DNmjHMf0qaN3xKlF7GYegTYF5Te5+UZKcgPP8DYsW74pil9I5lZvBjuuw/efNNvSdKPWFr8/wK+EpE3vPQ5mMkmZXnmGZgwAS6/3FpRRnLTsycsXEhZ9C6j9ojqlllEsoCBwC7geC/7M1WdkwDZyjC3zLWHKqxaZY4NjdRi8WJn+unZ029JUotIbpmjtvhVtVREHlPVI4Cv4yadEXeefx4GD4a2bU3pG6nFvn0wdKh7dj/7zG9p0oNYbPxTRWS4iE2iTlUKCpw/FPNia6QiderAs8+6xotRO8Ri478WuAnYJyJeXCbUm8SV1Czbsoyde3bSJqcNrbJbUScrMwert2kDM2bAQQf5LYlhVI9jjtm/vmOHDUyoKbE4aUvZSzz2i7E8Pst5lBCEVtmtaNu4LW1y2tA2p/xvm5w25bY1qtfIZ+lrTlERfPklnHIKHHaY39IYRs0ZNQqmToVZs6B+fb+lSV1imiIhIufiOncV17n7ZjyFqi1+0/83DMobxIbCDRQUFrBh5wYKigooKCxg5pqZFBQWsGPPjrD7Nq7feP/LoXFb2mRXfDkE0rkNc8mSWKNYJo5774W//AV+/NFa+0Z6cMop0Ly531KkPpUGWxeRx4GDgRe8rAuBJap6XZxlKyOeo3qK9haxsXBj2cuh7AVRWLD/heH9biraRKlWDO5cN6surbNbR/yaCOQHlvp1EtNUKSqCKVNcNC3DMDKPSKN6YlH8PwA9vChcgSGe36tqj7hIGoZkGc65r3Qfm4s3l3s5hHtBbNi5gQ2FG9hVsivscXIb5sb8NdGkfhOq2q/+ww+uhW+fwka68tVX8NRTLgJjVvJ9bCcN1RrO6bEY6ASs8NIdvbyMo05WnbJWe+82vaOWVVUK9xZW/HoIvDCK3Pr3Bd8zrXAaPxWHD1bWoE4D2kpj2qzbQZste2hLDm0GnEzbAadU6Jtold2K4sK6DBoEP/uZGwlhGOnIggXO19SqVTbBqzrE0uL/BOgPzMTZ+AfgomdtA1DVuBsSkqXFH0/27tvLpqJNFV4QG/77EQWfvc+GhvsoyIGCHNjQGPaGGaAkCC2zW9Jw4WUccNAWuhxSHNbcFMjLqR9r+GTDSC5UnSkzJ/QRnjwZRo+GlSvdhJUxYzI6Pm9NWvx3xkEeI4R6derRvkl72jdpX37DJY/Ain3lshTY1q0jBTOmlPuiWLqqiML6SyjovIKCwgK+Xue2bdu9LWydOfVyyr8UsiuamwLbWzRqkZQd2EZmIuKUfmkpvPeeC9soz092E1aKilyhFSv2R3HPYOUfjkpb/MlAJrT4I5KV5Zo3oYi4p97j2Wfht7+F6dMrDt3cVbKLjYUbo5qcAts2Fm5kn+4jlDpSh9Y5rSu+FCJ0YDes27C2r4RhVGDSJLjsMvjwQzj16jyn7EPp3BmWL0+0aElBTVr8hp906hT+YQ7xuzB4MFxzTXhfJg3rNqRjs450bNax0upKtZQtxVsqmpwKy3dmL9myhA07N1C4tzDscZo2aBr1BRG8LbdhbpU7sA0D4KKLXMt/8GCceScckfIzGGvxJzuTQz5fAbKz3XCG/HyKi6FhQ/cB4AeFewrZWLQxfCd2UfnhsZuKNqFhwjfXy6oXsR8i1OTUOrs19epYcGCjIkWdDiV71cKKG6zFX70Wv4g0AjqpapirasSVgG0yTIdVSYmLVNSjBzz2mD/i5dTPIad+Dnm5eZWWLSktYXPR5krnTMzfOJ8NOzewe9/usMdp0ahFpSanQLpx/cb2NZEBfPcdDNk+h4kNLmTI7nf2b8jONidVYYgl9OJZwFigPtBFRPoC9yRiNI/hkZ8ftnMqKwtOOAG6dvVBpmpQN6subRu3pW3jtpWWVVV27NkR1dxUUFjAdxu+o6CwgC27toQ9TsO6DWOaM5Hp/pxSnUMOgUGnN6Jdn2vhye9sVE8lxDKcczZwCvCx554ZEZmrqgnz/pLRpp4IlJbaxJVg9uzbE7EDO9wLo6S0pMIxIvlzitSJnV0v24czNYzYqYmpZ6+qbgv5XE7+joE0Zu5cN5LhhRecmceA+nXq06FpBzo07VBpWVVly64tlb4cYvHnFIvDvzY5bWjeqLkNh00AxcVw990wfDgMGOC3NMlLLIr/exG5BKgjIt2AG4Av4iuWEY2iIuejPDfXb0lSExGhRaMWtGjUgkNbHVpp+eK9xRXdc4S8MJZuWcqM1TMq9ecU1jVHyNdE6+zWNKjbIB6nnvbs3ev89jdtaoo/GrGYerKB0cAQXJD1D4A/qWp4RzRxwEw9FVH1bySPEZlgf07RHP4FthWXFIc9Tm7D3Ji/Jpo2aGod2EFs2wbNmvktRXJQbSdtyYApfsef/+wCUNxwgyn9dCDUn1O44bDB2zYXbw57nAZ1GkR8KYS+MFplt6JuVmZM31m+3LlwzuSXQLVt/CJyCPB7IC+4vKqeUpsCGtEpLYWvv3afsEZ6ICI0rt+Yxi0a07VF5UOzQv05hf2i2LmhbKTTnn17Ktbp+XOKpfO6bU7blPXnVFAAvXrB9dfDfff5LU3yEcur/xXgH8DTQMW5/EZCyMqCV1+FkhJr7WcqEf05hUFV2bZ7W6UuxOesn8OGnRsi+nPKrpcd9gUR7uuiRaMWSTMctk0bePhh56XWqEhMwzlV9cgEyROWTDf1fP65869/wAF+S2KkK7tLdlfagR34jeTPKUuyyjqwg+dNhIszYf6cEkOVTT0i0sJbfUdEfgO8AZRNpVTV8A7kjVqltBQuvxy6dHGxRg0jHjSo26BK/px+Kv4p4tdEYH3GqhlsLNrIzj07wx6naYOmEU1OoV8U1fXntGUL/PrX7j80dGiVd09bopl6ZuPG6weu9i1B2xSwKK4JICvLhU/cGf6/YxgJJ0uyaJXdilbZrejZOoxXwBBC/TmFe0ks+mkRn6/8PKo/p9Y5rSvMwg432a5NTpsyf045ObBwoflpCyUWU0/D0KGb4fLiSaabegwjU9hXuo9NRZtimoEdzZ9T84bN93dUN2pH2yZhXIp726sT3jRVqMnM3S+AfjHkGbXMa6/Bu+/CuHE2msfIDOpk1Snz53QY0b3CqCo79+ys9AUxb9N3TF2+gS2LukPuMmiyodxxGtZtGPOciXTx5xTNxt8O6AA0EpEj2G/yaQqYk5IEsGqV8zrYuLHfkhhG8iEiNGnQhCYNmnBwi4Ojlt2wATp1Ui6/upDf3Lyo7IuhzATlzZlYu2Mt36z/hoLCAvaW7q1Yp+fPKZKjv9AO7GT15xTR1CMiVwBXAkcB/2W/4t8OPKuqrydCQMhsU485YzOM2uH99+H442NrSKkqW3dtLWdWijQkdsPODTH7c4rmQjwe/pyqPXNXRIar6mu1Kk0VyTTFrwqLF0O3bn5LYhjpR2mpmwtTm2b9gD+nQPjScE7/Ai+PjUUbY/bn1Ca7DVf3u5oeravnjbHaNn6/lX4m8umnMGgQvPOOCyJtGEbtUFAAZ5/tZvReemntHbdRvUZ0zu1M59zOlZYt1VI2F20O67spOG/RZmeSOvOQM6ut+CORGU47UozDDoO//tWLI2oYRq3RqhW0bQuNGvknQ5Zk0TqnNa1zWtOLXpWWj4c/tWg2/vNV9RUR6aKqy2q95iqQaaYewzCM2iCSqSdaT8Lt3q+ZehLImDEwfbrfUhhGeqPq/PYvXeq3JP4QzdSzWUT+g4uz+3boRou5W/ts2waPPgp79sBxx/ktjWGkLwUFMHIkXHst/P3vfkuTeKIp/jNxk7SeA6p8aURkAvBzoEBVe3t5LYCXcC6elwMXqGr4KNkZSLNmrgWyz3ygGkZcadsWvvwSelbucSItiWjqUdU9qvolcKyqfoLz3TNbVT/x0pXxDHB6SN5twFRV7QZM9dIGLlaoqut0sglbhhF/evd2c2T2VAxbkPbEMlugrYjMAb4H5ovIbBHpXdlOqvopEOrBcxjwrLf+LHBOFWRNa266CU44wVr7hpFI5s9382WmTPFbksQSy3DO8cBNqvoRgIgM8vKOrUZ9bVV1nbe+HmgbqaCIjARGAnTq1KkaVaUWAwc6f/t1Ut8NiGGkDF27wpFHZt5Xdiwzd79V1T6V5UXYNw94N8jGv1VVc4O2b1HV5pUdx4ZzGoZhVJ3qDOcMsFRE/kdE8rzlj0B1B0FtEJH2nkDtgYJqHidt2LwZXnzRTDyGkVAmT4a8PGfkz8ujaMKLPPaYC22aCcSi+H8JtAZex43pb+XlVYe3gSu89SuAt6p5nLThuefg4otdsAjDMBLA5MluLOeKFW5ExYoVTPn1a1x/PXz4od/CJYZKTT3VPrDIC8Ag3ItiA3AX8CbwMtAJWIEbzllpCMd0NvWUlrphZcdWp8fEMIyqk5fnlH4QCsxu93OOWveOLyLFi2p750wG0lnxG4aRYLKyXEs/FBEoLaWkBOqmiRezmtj4jTiwa5cbyfNOejUwDCP5iTRKsFMnnn/eDe/cti2xIiWaShW/iFRwHhAuz6ga69e7hkemDSMzDN8ZMwayQyJjZWfDmDF07+6GdxYW+iNaoohlOOfXqtqvsrx4kq6mnsClT9M4z4aRvEyeDKNHw8qV7gtgzBjIz/dbqlqnyoFYROQY3CSt1iJyU9CmpoBNM6oB338PnTtba98wfCM/P6qiX7vWeck9//wEypRAopl66gONcS+HJkHLduC8+IuWnpSWwnnnuShAhmEkJ3/6E1x1Vfra+mMx9XRW1RVRC8WZdDP1fPmlcwx14ol+S2IYRjjWr3d2/q5d/ZakZlQ75i7QQETG41wpl5VX1VNqT7zMYuBAvyUwDCMa7drtX9+3L/18aMWi+F8B/gE8DZhjgRowbRpMnQq33272fcNIBW68EVatgtfSLA5hLOP4S1T1CVWdqaqzA0vcJUtDpk+HSZOgfn2/JTEMIxYOLPiaLlPGs0/quhm/kyf7LVKtEC3Yegtv9QacM7U3gN2B7bG4Wqgt0snGv3OntfYNIyUI+PQpKtqfl50N48enzNDPKrtsEJFlOBcW4UaZq6oeVLsiRiYdFP+WLdC8UgfUhmEkDUE+fb6hD3XYx2HMc2Oxly/3VbRYqXLnrqp2ia9ImcOcOXDMMfDGG3DGGX5LYxhGTKxcCcAe6nE673MMM3iDc8vyU5lKO3dF5Nww2duAuaqa8f70Y6FVK7jmGvPAaRgpRadOsGIF9dnL65xLT+bvz09xYhnVMwI4BvjISw/CBV7vIiL3qOpzcZItbejYER55xG8pDMOoEmPGlNn4j2WGy/N8+qQ6sYzqqQv0UNXhqjoc6Imz/R8N3BpP4dKBCRPgxx/9lsIwjCqTn+86cjt3BhGWHnA8p3ZZzJyeqdGxG41YFH9HVd0QlC7w8n4C9sZHrPRg+3YYNQoee8xvSQzDqBb5+a4jt7SUFt9/xpp97Vmzxm+hak4spp6PReRd3EQugOFeXg6wNV6CpQNNm8KiReZ90zDSgdxcmD8/Pf7PsSj+63DKPuCDfyLwmrpxoCfHS7BUR9U9IG3a+C2JYRi1hYj7b3/zDRxxhN/SVJ9KTT3qeFVVR3nLq5oK8Rp95o474NxznZ8PwzDShwcfhKOOSu2+u2j++D9X1eNFZAeuM7dsE+590DTu0qUwLVs6R0/p5tzJMDKd/Hxn9umSwjOdLNi6YRhGmlKjYOsicryIXOWttxKRFH7XxZfCQvjsM7+lMAwj3rz8MvzlL35LUT1iCbZ+F268/u1eVn1gUjyFSmUmTHABVr791m9JDMOIJx99BK+/DiUlOIdueXmQlZUSXjxjGdXzC+AI4GsAVV0rIk3iKlUKM2KEs+336eO3JIZhxJOxY6FRI8h6IcSL54oVLg1J68UzFlPPHm8UjwJ44/eNCGRnp2+AZsMw9pOT4xr4u+/4X3YUhajSoiIYPdofwWIgFsX/sog8CeSKyDXAFOCp+IqVepSUOIX/8cd+S2IYRqIoLoZDV37APdxZcWMSe/Gs1NSjqmNF5DRgO9AduFNVP4y7ZCnGihXw9dewdavfkhiGkSgaNYKRua8wYOsHFTcmsRfPWGz8eIrelH0UunaFhQvdp59hGJnD7Y92gJFfQlCgrmT34hltAlfoxK2yTdgErnKsXQtt20LdmF6jhmGkFfn57Ciuy5O3LGbE1r/TvHNTp/STtGMXokfgKhu5IyJzVDWFPVPED1X4xS9cWMX33/dbGsMw/GBp/wu5ZSu0eXY0l1/utzSVE2sbNfmn9/rIH/5grhkMI5Pp08eZeg85xG9JYsOMEzVEBIYP91sKwzD8JqD0A555k5loNv7gWLu5obF3VfX1uEmVIsyaBTNnuklbDRr4LY1hGH7z5JNu+e9/k9sKEK3Ff1bQ+ichaQUyXvG/8go8/TRcfrkpfsMw3Kz97t1h2zZo0cJvaSJj3jlrgCqsXu2CqRuGYSQbNfLOaVRk715nxzOlbxhGKKtXw+LFfksRGVP81WDxYqfwp071WxLDMJKNkhIXoevWW/2WJDI2qqcalJTAMcdAr15+S2IYRrJRty7885/Qs6ffkkQmJsUvIscCecHlVXVidSsVkeXADmAfUBLOBpXMHHoovPGG31IYhpGsnHmm3xJEJ5ZALM8BY4Hjgf7eUhuK+mRV7ZtqSv/f/4aNG/2WwjCMZGfpUrj6ati82W9JKhJLi/8ooKemwvCfOLNzJ1x4IVxwgfuUMwzDiERxMbz0knPX/rOf+S1NeWJR/POAdsC6WqxXgf+IiAJPqur40AIiMhIYCdApSdybNm7sJmw1auS3JIZhJDu9esG6dU5vJBuxKP5WwHwRmQnsDmSq6tk1qPd4VV0jIm2AD0XkB1X9NLiA9zIYD24cfw3qqlV69PBbAsMwUoWA0i8sdBG7koVYFP/dtV2pqq7xfgtE5A1gAPBp9L385cEH3TDOcePM/bJhGLEzahT85z8wd27yxOuIJQLXJ7VZoRezN0tVd3jrQ4B7arOOeFBQ4CZlmNI3DKMqnHQStGrlJn0mi2uXSl02iMhA4BGgB1AfqAMUVjcQi4gcBAQGQ9YFnlfVqKFqksVlQyp43TMMwwgQyWVDLO3XR4GLgFdwI3wuB6rtdVpVlwJ9qrt/otmzx8XT7dbNlL5hGNVDFT76CA44wM0D8puYLE6quhioo6r7VPVfwOnxFSt5mDTJ3ahvvvFbEsMwUpUdO2DYMHjoukWQl+eM/Xl5MHmyL/LE0uIvEpH6wDci8jfcsM4k6aKIP2eeCWPHugg7hmEY1aFpU/jw9+/T975LoHiLy1yxAkaOdOsJjs8bi42/M7ABZ98fBTQDHve+AhJCstj4DcMwqk1enlP2oXTuDMuXx6XKartlVtUVgADtVfV/VfWmRCp9v1CFm26COXP8lsQwjLRg5Uo+5zhOZhrbaVIuP9FUauoRkbNwvnrqA11EpC9wTw0ncCU9S5fCM884E88RR/gtTWLYu3cvq1evZteuXX6LYqQxDRs25MADD6RevXp+i5JYOnWiwYrdrKM9K+jMYcwry080sU7gGgB8DKCq34hIlzjKlBR07eq+yho29FuSxLF69WqaNGlCXl4eYkOYjDigqmzevJnVq1fTpUvaq5HyjBlD/5EjmV/Ukyw8E3t2NoyJOpo9LsTSSbtXVbeF5CWNC4V4sHOn+23SBDKpUbJr1y5atmxpSt+IGyJCy5YtM/OrMj8fxo8nq3MnSqjL2g79Yfz4hHfsQmwt/u9F5BKgjoh0A24AvoivWP5y7rnQrJkLpp5pmNI34k1GP2P5+ZCfzyknQp068FHidT4Qm+L/LTAa56DtBeAD4E/xFMpPVOGss8wDp2EY8eO3v4X69f2rPxZfPUU4xT86/uL4j4i7KYZhGPHi/PP9rT+i4heRt6PtmI6jen78ERYscC3+ZPGiZxhGerJ9O0ycCJdf7iZ4JZJoLf5jgFU4885XuLH8ac0TT7i+lpUroWVLv6UxMoE333yT9957j+3btzNixAiGDBnit0hGgli40FkXWraEiy9ObN3R2rXtgDuA3sA44DRgk6p+UtuumpOF+++Hzz4zpe83jcOELDr22GMTLkedOnXo27dv2bLcm1157LHHsnXrVh5//PGYj/Xkk0/Srl07+vTpQ9euXZk4cSIA55xzDk899RT/+Mc/eOmll6ot6/vvv0/37t05+OCDuffee2Mus2rVKk4++WR69uxJr169GDduXLVlMKpG//7OR3+ilT7gxtVWtgANgCuBjcD1sexTm8uRRx6p8aa0NO5VJD3z58/3WwRVVc3JyUlofaWlpbpv374qybFs2TLt1atXzHVcd911+sQTT6iq6ldffaUtW7Yst/2mm27S2bNnx3y8YEpKSvSggw7SJUuW6O7du/Xwww/X77//PqYya9euLat3+/bt2q1btwr7xoNkedbSHWCWhtGpUS3ZItJARM4FJgHXAQ+z35d+2rBunZuhO32635IYkQh8BSxfvpwePXpwzTXX0KtXL4YMGUJxcTEAkyZNYsCAAfTt25drr72Wffv2Aa5VfeSRR9KrVy/Gjx9fdpzu3btz+eWX07t3b1atWlUlWW677TaWLFlC3759ueWWWyrd57vvvqN79+4AdOnShfrekA5V5dZbb+WMM86gX79+sV+QIGbOnMnBBx/MQQcdRP369bnooot46623YirTvn37snqbNGlCjx49WLNmTbXkMKrHvffCNdckts5onbsTcWaefwP/q6rzEiZVgtmwwc3Qbd/eb0mMWFi0aBEvvPACTz31FBdccAGvvfYaRx55JC+99BLTp0+nXr16/OY3v2Hy5MlcfvnlTJgwgRYtWlBcXEz//v0ZPnx42XGeffZZBg4cGLae4uJi+vbtCzhl/cYb+9s89957L/PmzeObGP11z507l+7du6OqPProo4zxZms+8sgjTJkyhW3btrF48WJ+9atfAXDCCSewY8eOCscZO3Ysp556arm8NWvW0LFjx7L0gQceyFdffVXlMsuXL2fOnDkcffTRMZ2TUTsUFrqO3tLSxA0qida5eylQCPwOuCFo0oUAqtWMwJWM9O0LM2f6LUVyceP7N/LN+m9q9Zh92/XlodMfqvFxunTpUqaQjzzySJYvX87WrVuZPXs2/fv3B5zSbtOmDQAPP/xwmdJetWoVixYtol27dnTu3Dmi0gdo1KhRzIo9GqtWrWLHjh0MHTqUNWvWcPjhh3P33XcDcMMNN3DDDTdU2Oezzz6rcb1VYefOnQwfPpyHHnqIpokeYpLh3HNP4oM8RVT8qpoRAxpnzYLevTPLJ0+q0yAocGmdOnUoLi5GVbniiiv461//Wq7sxx9/zJQpU5gxYwbZ2dkMGjSozF1ATk5OQuSdO3cuJ554ItOmTWPLli307t2bGTNmRO2wrkqLv0OHDuVMVatXr6ZDhw4xl9m7dy/Dhw8nPz+fc889t1rnaFSfgNLfsAFycxMTlzejQ4cXFsKQIXD22c4Tp7Gf2miZJ5LBgwczbNgwRo0aRZs2bfjpp5/YsWMH27Zto3nz5mRnZ/PDDz/w5Zdf1kp9TZo0qaCYBw8ezMSJEyso3e+++44jPBevzZs355JLLuG9996Lqvir0uLv378/ixYtYtmyZXTo0IEXX3yR559/PqYyqsqIESPo0aMHN910U8x1GrXLvHnQrx9MmACXXhr/+jKiVR+J7Gx47TX4/e/9lsQIpqioiAMPPLBseeCBByrdp2fPnvz5z39myJAhHH744Zx22mmsW7eO008/nZKSEnr06MFtt90W1bRTFVq2bMlxxx1H7969ueWWWygtLWXx4sW0aNGiQtm5c+eWKX6As846i3//+9+1IgdA3bp1efTRR/nZz35Gjx49uOCCC+jVqxcAQ4cOZe3atRHLTJ8+neeee45p06aVDVutTdmM2OjZE0aPhlp6PCul0ghcyYBF4EoMCxYsoEePHn6LkZLMmzePCRMmxPSSMuxZSxTVjsCVrjz3nHODXVLityRGOtC7d29T+kaNmTvX6aZ4k7GK/4sv4L33nGtUwzCMZOCxx+B3v4N4hyvIWMX/xBMwdWrih1EZhmFE4q67YNGi+I8yzLhRPaWlsGWL88djPvcNw0gmEjWJNONa/G++6WIbz5njtySGYRgVWb/euYaP5+CqjFP8vXvDr34Fhx/utySGYRgVadkS1qxxlol4kXGmnkMOgb//3W8pDMMwwlOvHsyeHd/+x4xq8T/wAHgu1Q3DMJIWERf/e+XK+Bw/YxT/0qVw220Q4q3WMAwjKbntNucuvqio9o+dMaaegw5yyj/MjHrDMIyk44ILoGvX+LhqzgjFH/BzfeCBfktiGPHBYvemH0ce6ZZ4kBGmnvx8uPZav6UwYmX16tUMGzaMbt26cdBBB3H99deze/fuWjn23XffzdixY2MuHy62blXj/6Zb7N5x48bRu3dvevXqxUMPPVRu24MPPkivXr3o3bs3F198cZkLbCPJCBePMdmWmsTcLS1VveMO1TFjqn2IjCEZ4qCWlpZq//79dcKECarqYsX+8pe/1BtuuKFWjn/XXXfp/fffX6HOcDF3VaseWzcc6RS7d+7cudqrVy8tLCzUvXv36uDBg3XRokWqqrp69WrNy8vToqIiVVU9//zz9V//+lfYupLhWcsEqE7M3XRAxDlju+MOvyUxYmHatGk0bNiQq666CnCt5QcffJCJEycyb948evfuXVZ27NixZZGsIHxsXYAxY8ZwyCGHcPzxx7Nw4UIgfMzdcPuHi60biP8LMHHiRA4//HD69OnDZZddVqVzTcXYvQsWLODoo48mOzubunXrctJJJ/H666+XbS8pKaG4uJiSkhKKioo44IADqiWLEWfCvQ2Sbalui3/1atWZM6u1a0YS2go76STVQINtzx6Xfu45ly4sdOkXX3TprVtd+rXXXHrjRpd++22XXrcuNhnGjRunN954Y4X8vn376pw5c8q1ju+//3696667ytKbN29WVdWioiLt1auXbtq0SWfNmqW9e/fWwsJC3bZtm3bt2lXvv/9+XbZsmYqIzpgxI+r+4VrkgRb8vHnztFu3brpx48Zy+4eSlZWlffr00T59+ug555xT7jhVbfHn5ubqmjVrtLS0VO+8886yL6Nx48Zpv3799Nprry37IlBVPf7448vqDl4+/PDDCsd+5ZVXdMSIEWXpiRMn6nXXXVeuzPz587Vbt266adMmLSws1IEDB+r1119ftv2hhx7SnJwcbdWqlV5yySURz8Na/ImBCC3+tO7cXbUKRo6EKVPAC79qpDHhYut++eWX/OIXvyA7OxuAs88+u6x8aMzdSLF5IzFt2jTOP/98WrVqBRA2CAukV+zeHj16cOuttzJkyBBycnLo27cvdTwXt1u2bOGtt95i2bJl5Obmcv755zNp0iQuTURIKaNKpLXiHzgQvv3WPHBWl48/3r9er175dHZ2+XSzZuXTrVqVT0fRn+Xo2bMnr776arm87du3s379elq2bElpaWlZfnDHYbTYupEIjrlbnf0TTTLE7gUYMWIEI0aMAOCOO+7gQG+43JQpU+jSpQutW7cG4Nxzz+WLL74wxZ+EZISN30gdBg8eTFFRUdlolX379nHzzTdz/fXX065dOwoKCti8eTO7d+/m3XffLdsvUmzdE088kTfffJPi4mJ27NjBO++8E7beSPuHi60b4JRTTuGVV15h8+bNAPz0009VPt9IsXvXrFlToWyk2L3R+Oyzz/jmm28qLKFKH8rH5d2zZw8vvvhiuS+kAAUFBQCsXLmS119/nUsuuQSATp068eWXX1JUVISqMnXqVIuylaT4ovhF5HQRWSgii0XktrhUMnky5OW5Afx5eS5tJD0iwhtvvMGrr75Kt27daNmyJVlZWYwePZp69epx5513MmDAAE477TQOPfTQsv0ixdbt168fF154IX369OGMM86gf//+YeuNtH9obN1gevXqxejRoznppJPo06dPtYKVp1rsXoDhw4fTs2dPzjrrLB577DFyc3MBOProoznvvPPo168fhx12GKWlpYwcObLW5Ms44qnDwhn+47kAdYAlwEFAfeBboGe0farcuTtpkmp2tqpzd+GW7GyXb0QkGTvcpk+frp06dar2EMVUY+7cuTpq1Ci/xYg7yfisJRW1pMOI0Lmb8GDrInIMcLeq/sxL3+69gP4aaZ8qB1vPy4MVKyrmd+5sXtqiYAGwjURhz1ol1JIOS6Zg6x2AVUHp1V5eOURkpIjMEpFZGzdurFoNkVzaxcvVnWEYRm0SZx2WtJ27qjpeVY9S1aMCowRiplOnquUbhmEkE3HWYX4o/jVAx6D0gV5e7TFmjBtvGEx2tss3DMNIduKsw/xQ/P8FuolIFxGpD1wEvF2rNeTnw/jxzh4m4n7Hj3f5RlQS3edjZB72jMVAnHVYwjt3AURkKPAQboTPBFWN+hqrcueuUS2WLVtGkyZNaNmyJWITIIw4oKps3ryZHTt20KVLF7/FSXside76MnNXVf8NxDGGvFEdDjzwQFavXk2VO9MNowo0bNiwbLav4Q9p7bLBqBr16tWzVphhZABJO6rHMAzDiA+m+A3DMDIMU/yGYRgZhi+jeqqKiGwEwsxfjolWwKZaFCcZSfdztPNLfdL9HJP1/DqraoUZsCmh+GuCiMwKN5wpnUj3c7TzS33S/RxT7fzM1GMYhpFhmOI3DMPIMDJB8Y/3W4AEkO7naOeX+qT7OabU+aW9jd8wDMMoTya0+A3DMIwgTPEbhmFkGGmt+BMS1D3OiEhHEflIROaLyPci8jsvv4WIfCgii7zf5l6+iMjD3jl/JyL9/D2D2BCROiIyR0Te9dJdROQr7zxe8lx4IyINvPRib3uer4LHiIjkisirIvKDiCwQkWPS6R6KyCjv+ZwnIi+ISMNUv4ciMkFECkRkXlBele+ZiFzhlV8kIlf4cS6hpK3iF5E6wGPAGUBP4GIR6emvVNWiBLhZVXsCA4HrvPO4DZiqqt2AqV4a3Pl285aRwBOJF7la/A5YEJS+D3hQVQ8GtgAjvPwRwBYv/0GvXCowDnhfVQ8F+uDONS3uoYh0AG4AjlLV3jh36xeR+vfwGeD0kLwq3TMRaQHcBRwNDADuCrwsfCVcBPZ0WIBjgA+C0rcDt/stVy2c11vAacBCoL2X1x5Y6K0/CVwcVL6sXLIuuChsU4FTgHcBwc2CrBt6L4EPgGO89bpeOfH7HCo5v2bAslA50+Uesj+OdgvvnrwL/Cwd7iGQB8yr7j0DLgaeDMovV86vJW1b/MQY1D2V8D6JjwC+Atqq6jpv03qgrbeeiuf9EPAHoNRLtwS2qmqJlw4+h7Lz87Zv88onM12AjcC/PHPW0yKSQ5rcQ1VdA4wFVgLrcPdkNul1DwNU9Z4l5b1MZ8WfVohIY+A14EZV3R68TV1TIiXH5YrIz4ECVZ3ttyxxpC7QD3hCVY8ACtlvIgBS/h42B4bhXnAHADlUNJGkHal8z9JZ8cc/qHuCEJF6OKU/WVVf97I3iEh7b3t7oMDLT7XzPg44W0SWAy/izD3jgFwRCQQKCj6HsvPztjcDNidS4GqwGlitql956VdxL4J0uYenAstUdaOq7gVex93XdLqHAap6z5LyXqaz4o9/UPcEICIC/BNYoKoPBG16GwiMELgCZ/sP5F/ujTIYCGwL+jRNOlT1dlU9UFXzcPdomqrmAx8B53nFQs8vcN7neeWTutWlquuBVSLS3csaDMwnTe4hzsQzUESyvec1cH5pcw+DqOo9+wAYIiLNvS+jIV6ev/jdyRDPBRgK/AgsAUb7LU81z+F43Ofkd8A33jIUZxOdCiwCpgAtvPKCG820BJiLG2nh+3nEeK6DgHe99YOAmcBi4BWggZff0Esv9rYf5LfcMZ5bX2CWdx/fBJqn0z0E/hf4AZgHPAc0SPV7CLyA67PYi/tqG1Gdewb80jvXxcBVfp+XqprLBsMwjEwjnU09hmEYRhhM8RuGYWQYpvgNwzAyDFP8hmEYGYYpfsMwjAzDFL9RK4jIzmrud6WIPFrLsjQQkSki8o2IXFibxw6pJy/Yc2OiiPWaicjH4rzTfisi04PmEYQre4CIvBrDMe+oqrxG8mGK30hHjgBQ1b6q+lLwBs9rayaRr6p9gGeB+yMVUtW1qnpepO1BmOJPA0zxG7WKiAzyWpoB3/OTvdmciEh/EfnCa4HOFJEm3m4HiMj7nr/yvwUda4iIzBCRr0XkFc9fESJyr7j4BN+JyNiQ+tsAk4D+Xou/q4gsF5H7RORr4HwRuVhE5orzHX9f0L47ReR+cX7lp4jIAO9clorI2VW4BteIyH+983xNRLK9/GdE5Lygcjtr+5pF4VPgYG9m6f3euc8NfBEFf714XxSvhx5fRO4FGnnXdXKs18NIQvyeQWZLeizATu93EM7b4oG4hsUM3Ozj+sBSoL9XrinOedmVXn4z3IzOFTjfJq1wyirHK38rcCdu5uRC9seLzg0jyyC8GcBeejnwB2/9AJyLgdZe/dOAc7xtCpzhrb8B/Aeoh/Of/02YevIIctkblN8yaP3PwG+99WeA8+J1zcLI8THeDFLgFuAlYDjwIc5nflvvWrQPPpdoxw/IbEtqLwEHSoZRm8xU1dUAIvINTqlsA9ap6n8B1PMw6jVsp6rqNi89H+gM5OIC6Ez3ytTHKcRtwC7gn+Kidb0bo0wBk09/4GNV3ejVNxk4EedGYQ/wvlduLrBbVfeKyFzvHGKlt4j82TuHxsTmm6U2rtmqCkeFySJSjHv5/Ra4CXhBVffhHI59grsm34XsF+vxjRTEFL8RD3YHre+j8ucsXHkBPlTVi0MLi8gAnCOw84DrcR49K6MwhjJ71WvW4mID7AZQ1VLZ72UyFp7BfUV8KyJX4lr04KKpZQGISBbuZRagNq5ZOPJVdVYg4b00YqGq8hgphNn4jUSxEGgvIv0BRKRJJcr0S+A4ETnYK58jIod4dv5mqvpvYBTODFMVZgIniUgrr6P3YuCTqp5MJTQB1olzp50flL8cONJbPxtnRopGVa9ZLHwGXCguxnFr3NfOzCrsv9c7LyOFsbe4kRBUdY/XkfiIiDQCinF+3COV3+i1ll8QkQZe9h+BHcBbItIQ91VwUxXlWCcit+FcBgvwnqq+Vclu0eguIquD0qOA/8FFSdvo/QY6ZJ/yZP8WZ1KK+hVS1WsWI2/gwiB+i+vT+IOqrpfYA56PB74Tka/Vuc82UhDzzmkYhpFhmKnHMAwjwzDFbxiGkWGY4jcMw8gwTPEbhmFkGKb4DcMwMgxT/IZhGBmGKX7DMIwM4/8BP+FubwADviAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "processTrajectories(d2hs,trials)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This tells us that \n",
    "\n",
    "* 1 less than 2% of the variation in the measured data can be explained by the linear model,\n",
    "\n",
    "* 2 more than **98%** of the  variation can be explained by the quadratic model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Matplotlib.pyplot: Formatting text: LaTeX\n",
    "```python\n",
    "  plt.plot(distances, altitudes, 'g', label = \"Linear Fit, $R^2$= {:.2f}\".format(r1line))\n",
    "```\n",
    "Matplotlib has great support for LaTeX. All we need to do is to use dollar signs encapsulate LaTeX in any text (legend, title, label, etc.).\n",
    "\n",
    "* For example, `$R^2$`. $R^2$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  4 Python Script\n",
    "\n",
    "```\n",
    "python trajector.py\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting ./trajector.py\n"
     ]
    }
   ],
   "source": [
    "%%file ./trajector.py\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def rSquared(measured, predicted):\n",
    "    \"\"\"Assumes measured a one-dimensional array of measured values\n",
    "               predicted a one-dimensional array of predicted values\n",
    "       Returns coefficient of determination\"\"\"\n",
    "    # RSS: residual sum of squares\n",
    "    estimateError = ((predicted - measured)**2).sum()\n",
    "    \n",
    "    # TSS: total sum of squares\n",
    "    meanOfMeasured = sum(measured)/float(len(measured))\n",
    "    variability = ((measured - meanOfMeasured)**2).sum()\n",
    "    \n",
    "    return 1 - estimateError/variability\n",
    "\n",
    "def getTrajectoryData(fileName):\n",
    "    dataFile = open(fileName, 'r')\n",
    "    d2hs=[]\n",
    "    discardHeader = dataFile.readline()\n",
    "    for line in dataFile:\n",
    "        dhs = line.split()\n",
    "        d2h={'d':None,'h':[]}\n",
    "        # the distance in first column\n",
    "        d2h['d']=float(dhs[0])\n",
    "        trials=len(dhs)-1\n",
    "        for i in range(trials):\n",
    "            d2h['h'].append(float(dhs[i+1]))\n",
    "        d2hs.append(d2h)\n",
    "    dataFile.close()\n",
    "    return  d2hs,trials\n",
    "\n",
    "def processTrajectories(d2hs,trials):\n",
    "    numTrials = trials\n",
    "    distances=[]\n",
    "    meanHeights=[]\n",
    "    \n",
    "    for item in d2hs:\n",
    "        distances.append(item['d'])\n",
    "        meanHeights.append(np.mean(item['h']))\n",
    "\n",
    "    plt.title('Trajectory of Projectile (Mean Height of {} Trials'.format(trials))   \n",
    "    plt.xlabel('Inches from Launch Point')\n",
    "    plt.ylabel('Mean Height of the projectile')\n",
    "   \n",
    "    # the experimental data\n",
    "    plt.plot( distances, meanHeights, 'ro')\n",
    "   \n",
    "    # Linear Fit\n",
    "    a,b = np.polyfit(distances, meanHeights, 1)\n",
    "    altitudes = a*np.array(distances) + b\n",
    "    \n",
    "    r1line=rSquared(meanHeights, altitudes)\n",
    "    plt.plot(distances, altitudes, 'g', label = \"Linear Fit, $R^2$= {:.2f}\".format(r1line))\n",
    "    \n",
    "    # Linear Fit\n",
    "    a,b,c = np.polyfit(distances, meanHeights, 2)\n",
    " \n",
    "    altitudes = a*(np.array(distances)**2) +  b*np.array(distances) + c\n",
    "             \n",
    "    r2quad=rSquared(meanHeights, altitudes)          \n",
    "    plt.plot(distances, altitudes, 'b:',label =\"Quadratic Fit, $R^2$= {:.2f}\".format(r2quad))\n",
    "    plt.legend()\n",
    "    plt.show()\n",
    "\n",
    "if __name__ == \"__main__\":  \n",
    "    fileName='./data/projectileData.txt'\n",
    "    d2hs,trials = getTrajectoryData(fileName)\n",
    "    processTrajectories(d2hs,trials)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Further Reading : When Theory Is Missing\n",
    "\n",
    "In an ideal world, we would run **a controlled experiment** (e.g., hang weights from a spring), study the results, and retrospectively formulate a model consistent\n",
    "with those results.\n",
    "\n",
    "Unfortunately, in many cases it is **impossible to run even one controlled experiment**\n",
    "\n",
    "In such situations, one can **simulate a set of experiments** by dividing the existing data into **a training set** and **a holdout set**.\n",
    "\n",
    "### How does one choose the training set?\n",
    "\n",
    "1. One way to do this is to **randomly choose the samples** for the training set.\n",
    "\n",
    "2.  A related but slightly different way to check a model is to train on **many randomly selected subsets** of the original data, and see how similar the models are to one another. \n",
    "\n",
    "If they are quite similar, than we can feel pretty good. This  approach is known as **cross validation**. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.7"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": false,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {
    "height": "calc(100% - 180px)",
    "left": "10px",
    "top": "150px",
    "width": "165px"
   },
   "toc_section_display": true,
   "toc_window_display": true
  },
  "widgets": {
   "state": {},
   "version": "1.1.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
